Positioning accuracy improvement in high ‐ speed GPS receivers using sequential extended Kalman filter

It is well known that in high ‐ speed movements, the positioning accuracy of global positioning system (GPS) receivers decreases drastically. The models presented so far to describe high ‐ speed motion do not represent the state of the system precisely, so positioning accuracy with the methods based on these models is not appropriate. Here, a comprehensive method is proposed to solve the accuracy issue of the single ‐ frequency GPS receiver at high ‐ speed motions without increasing the computational complexity. Suitable modelling of the GPS receivers at high ‐ speed motion, using the sequential extended Kalman filter, correct determination of the process noise covariance matrix and accurately estimating the variance of the observations are the basics of the proposed approach. Simulations with different data and motion scenarios (at speeds from 100 to 7300 m/s) show that the proposed method, while not increasing the computational cost, improves the accuracy of positioning more than 70% when compared to the conventional methods.


| INTRODUCTION
Nowadays, with the increasing use of global positioning system (GPS) in various fields, many efforts are being made to improve its performance. One of the main parts of a GPS receiver is the navigation section. In this section, navigation equations are formed using the pseudo-range and the carrier phase data from at least four satellites; then, by solving these equations, the user position is obtained [1,2].
By increasing the speed of the receiver, the accuracy of the positioning decreases. This is due to the jump in the output data caused by sudden acceleration, increasing the noise of the observations and decreasing the data quality. However, precise positioning is essential in high-speed movements such as aircraft or satellite motion [3]. Accurate positioning and parameter estimation in such dynamic stress is very challenging [4]. The methods presented so far, including the differential positioning, use of auxiliary tools and combination with other navigation systems, are either not sufficiently accurate at such speeds or have high complexity and implementation costs [5,6].
Two main methods of solving navigation equations are the least squares (LS) and the Kalman filter (KF). In the LS method, the position is obtained by minimizing the sum of the squares of the residuals of every single equation. The advantage of the LS method is its simplicity, but this method has less accuracy, and it also depends on the starting point. In addition, the performance decreases at high speeds or high noise levels. There have been several attempts to improve the performance of this method, including partitioned LS, recursive LS (RLS), constrained LS, weighted LS, recursive weighted LS (RWLS), non-linear LS and so on [5,7]. Although these methods improve the accuracy of the conventional LS, it is still insufficient to analyse a high-speed motion because the system model and movement parameters are not applied in these methods. The more complex methods may have better accuracy; however, they result in longer run times and tremendous computational cost.
Applying the KF algorithm is another way to solve the navigation equations, as mentioned before. In terms of accuracy, this method is better than the LS method. Also, it can overcome the influence of observation noise error better than the LS method if the system is modelled correctly. Another advantage of the KF over the LS method is that it identifies a trend of the receiver movement parameters [8]. As another point, the KF algorithm requires a dynamic model that can accurately describe the motion of the target. KF depends on the quality of this model to provide a priori knowledge of the system dynamics [9]. Some of the generalizations of the KF algorithm are extended KF (EKF), unscented KF (UKF), robust blind KF, adaptive fading KF, fuzzy adaptive cubature KF, Tobit KF and so on. [10][11][12][13][14][15]. All of these methods, mainly tested at not too high speeds, employ a simple model for accelerated movements, which results in a reduction in their performance at high-dynamic conditions.
Compared to the LS method, the KF is employed in positioning under a smaller user's dynamic range, that is, conditions in which acceleration's derivative is not more than 2g (g is the gravitational acceleration). When the user is under high-dynamic conditions, adopting the KF position and velocity results leads to a decrease in the performance. This degradation is caused by the estimation inaccuracy of the highest order's state, which arouses the model mismatch at the moment of acceleration's derivative appearing [16]. The addition of auxiliary tools such as a fuzzy system to improve the performance, although it may increase the accuracy, it also increases the computational cost and complexity.
Over the past years, various models of target motion have been developed [17]. There have also been several interesting studies on positioning using KF. Nevertheless, they are scattered throughout the literature, and there is no comprehensive study concerning how to choose a suitable model and apply effective positioning techniques on KF to increase accuracy without increasing computational complexity.
Here, by an appropriate modelling of receiver motion, sequential EKF, correct determination of process noise covariance matrix and proper selection of variance estimation method, a method is proposed to increase the accuracy of single-frequency GPS receiver positioning at high speeds. As it is shown in the paper, six data sets with different movement scenarios and speeds ranging from 100 to 7300 m/s are used to test and analyze the performance of the proposed method. It is also compared with LS, RWLS, and conventional EKF methods. The result of the simulations shows that the proposed method, while decreasing the computational cost, improves the positioning accuracy by more than 70% when compared to the conventional methods.
This study is organized as follows: Section 2 introduces the navigation equations. Section 3 describes the solution of the navigation equations using the EKF algorithm. In Section 4, a suitable model for high-dynamic movement analysis is presented and based on it, the transition and process noise covariance matrices are obtained. Section 5 discusses the selection of the proper method for estimating the variance of observations. Section 6 describes the sequential EKF to reduce computational complexity. Section 7 presents the results of implementing the EKF algorithm based on the model proposed in Section 4 and compares it with LS, RWLS, and conventional EKF methods. Finally, Section 8 presents the conclusions.

| NAVIGATION EQUATIONS
Positioning in the GPS is performed using the information from at least four satellites. Having the pseudo-range between a satellite and a receiver, the pseudo-range equation can be formed as follows [18]: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where ρ j is the measured pseudo-range from satellite j to the receiver, (X j , Y j , Z j ) is the j-th satellite position, (x, y, z) is the receiver position, c is the speed of light and t r is the receiver clock bias. By linearizing Equation (1) around the initial point (x 0 , y 0 , z 0 , (ct r ) 0 ), the following equation is obtained: So, Equation (1) can be rewritten as below: Δy ¼ y − y 0 ð9Þ By forming Equation (3) for at least four satellites, the following matrix equation is established:

| EXTENDED KALMAN FILTER ALGORITHM
The KF was introduced by Kalman in 1960 and is a powerful estimator under random noises distributions and mathematical model conditions [19]. It is a recursive filter that estimates the state of a dynamic system from a series of incomplete and corrupted evaluations [20]. It does not need to save previous information at each step and respond quickly; so, it is an excellent choice for real-time applications. An EKF is used as the estimator in this work. A general non-linear time-varying system can be described by the following equations: The state equation: The measurement equation: where X k is the system state vector at time t k , f k is the transition matrix at t k , that is, a non-linear function that defines the system's dynamics, Z k is the observation vector at t k and h k is the measurement matrix at t k , that is, a matrix that shows the ideal (noiseless) connection between measurements and the state vector. W k and V k represent the process noise and the measurement noise at t k , respectively, which both assumed to be additive, white, zero mean and uncorrelated with each other.
To obtain the best estimate of X k+1 , the EKF filtering equations are stated below: Estimation step: where F k and H k are the Jacobian matrices of f k (.) and h k (.) at t k , respectively, K k+1 is the Kalman gain, P k|k is the estimation is the process noise covariance matrix at t k .

| HIGH-DYNAMIC GPS RECEIVER MODEL
For a high-dynamic receiver, the conventional positionvelocity-acceleration (PVA) model, which is based on the Wiener process acceleration model, is as below [21]: In this model, the acceleration is modelled as a Wiener process, that is, the integral of the white noise. The Wiener process relates to the zero-mean white noise n(t) as follows [22]: where w(t) represents the Wiener process. Using this model, the state vector and the transition matrix for a high-speed receiver are defined as below: where X k is the state vector at t k , v x is the receiver speed along the x-axis, a x is the receiver acceleration along the x-axis and T is the sampling period. Equation (23) is the same for y and z components. Also, the measurement matrix and the process noise covariance matrix for a high-speed receiver are defined as follows [17]: where ρ i , i = 1, 2, …, N is the measured pseudo-range from satellite i to the receiver, S p is the spectral amplitude for the position random process and S f and S g represent the spectral amplitudes for the clock model. The measurement noise covariance matrix is given below: where σ 2 x n is measurement noise variance of the n-th observation and σ xixj is measurement noise variance between i-th and j-th observation. When the observations are independent, all the off-diagonal terms will be zero; therefore, the covariance matrix will be diagonal: In stochastic modelling, an unknown time-invariant quantity is represented by a random variable, while an unknown time-varying quantity is modelled by a random process. The simplest class of random processes is white noise. The second simplest class is either the Wiener process (the process with independent increments) or the Gauss-Markov process, which contain white noise and Wiener processes as special cases. White noise is isolated in the time because its value at one time is uncoupled from any other time, while a Gauss-Markov process is local in the time since its value at one time depends on the values at its nearest neighbours. Therefore, it is natural to consider a Gauss-Markov process model if white noise and Wiener process models are not good enough [17].
A stationary process that satisfies the requirements of the Gaussian process and has an exponential autocorrelation function is called a Gauss-Markov process. The autocorrelation and power spectral density (PSD) functions for a stationary firstorder Gauss-Markov process are as follows: where σ 2 is the variance of the process and 1/β is its time The process is non-deterministic, so a typical sample time function does not show any deterministic structure and looks like typical 'noise'. For two reasons, the Gauss-Markov process is important in practical applications: (1) it seems to fit a large number of physical processes with reasonable accuracy and (2) it has a relatively simple mathematical description [21].
Accordingly, the PVA model for a high-dynamic receiver based on the Gauss-Markov model is as follows: Although the acceleration is considered stochastic, this model is an approximation at best. There is no exact way to incorporate deterministic accelerations into any of the randomprocess dynamics appropriate for the KF model. Using this PVA model, the KF gives a better estimation than the PVA model presented in Figure 1, but it can still be inadequate.
To define the state transition and process noise covariance matrices in this model, first, the random process differential equation must be written: where w is the white noise and G is a matrix that gives the appropriate linear connections for the process. For the model drawn in Figure 2: x v x a x where w(t) is the unity white noise (e.g. a white noise whose PSD is unity). Considering the Laplace transform of Equation (33), ϕ k, which is the discretization of the continuous-time state transition matrix can be yielded: where s is the Laplace transform frequency parameter. Therefore: 1 Finally, it is obtained: The Equation (36) is the same for y and z components, and the clock model is the same as Equation (24). The process noise covariance matrix is also obtained as Equation (38). In this case, E[u(η)u(ζ)] is a Dirac delta function. Thus: So, by putting ϕ (from Equation (37)) and G (from Equation (34)) in Equation (39), Q xyz is obtained (see Equation (40). In the final form, Q is a block covariance matrix similar to Equation (26).

| VARIANCE ESTIMATION METHOD
The loss of quality of the measurements distorts them; As a result, the positioning accuracy is reduced. The loss of quality is mainly due to the random noise of observations, physical correlations, and unmodelled systematic effects [23,24]. These effects are measured by variance estimation methods and applied to the calculations by the R matrix in the KF algorithm. The two main factors used to determine the variance of observations are the elevation angle and the signal-to-noise ratio of the observations. So far, several relations based on these two factors have been proposed to estimate the variance, such as exponential function, sinusoidal function, tangent function, sigma ε method, sigma Δ method, C/N0-based function combined with noise bandwidth of the carrier tracking loop and so on. We have collected these models from different references and tested them in [5,25] and based on test results; it is found that the best variance estimation model is the exponential model given by Euler and Goad in [26]. So, we used it for our proposed method. This model is as follows: where σ 2 x k ðnÞ is the variance of the n-th observation at epoch t k , φ k (n) is the elevation angle of the n-th satellite at epoch t k , φ 0 is the reference elevation angle and the parameters a and b are obtained using the estimated variance components.

| COMPUTATIONAL COMPLEXITY
One of the significant criteria for evaluating an algorithm is its computational complexity. The importance of the issue manifests when the algorithm is used for high-dynamic receiver's positioning. The computational complexity of an algorithm is the number of resources required to run it. The number of arithmetic operations is one of the resources that is commonly used. It is evaluated using floating-point operations per iteration (flops), which is defined as addition, subtraction, multiplication, or division operations between two floating numbers [27]. For example, in the summation of two (n � m) dimension matrices, nm flops are generated in the computer. The matrix inversion in calculating the KF gain (Equation 18) causes a considerable computational cost, so eliminating this step improves the speed of the KF algorithm. For this reason, we propose the technique of performing EKF calculations sequentially (SEKF). If the observation matrix is a vector and the observations are not correlated with each other, that is, the measurement noise covariance matrix is diagonal, the calculations related to scalar observations can be performed separately and sequentially. This technique has two advantages: (1) The number of arithmetic calculations required for processing measurements sequentially is significantly less than the corresponding number for processing vector measurements. (2) Avoiding the matrix inversion by making the expression H k+1 P k +1|k H k+1 T + R k+1 a scalar, causes a significant reduction in the number of mathematical operations [28]. The expressions required to process EKF sequentially are as follows: where i = 1, 2, …, l and l is the number of observations, H k+1,i is the i-th row of H k+1 , R k+1,i is the i-th diagonal element of R k+1 and Z k+1,i is the i-th observation. The initial and final values are as follows:X kþ1jk;0 ¼X kþ1jk ð46Þ

| SIMULATION RESULTS
The Rohde and Schwarz GNSS simulator is used to generate raw GPS data such as pseudo-range, Doppler shift, integrated carrier phase and satellite ephemeris. The GNSS simulator in the vector signal generator supports all possible scenarios, from simple setups with static satellites to flexible scenarios generated in real-time with up to 24 dynamic GPS satellites. It can simulate realistic transmission conditions by modelling various atmospheric effects and using the multi-path signal generator. The simulator is used to model the effects that impact GPS receiver performance, such as multi-path reflections, atmospheric conditions, interference signals, and antenna characteristics. Six different motion scenarios have been used to test the proposed method. The first scenario involves a moving aircraft at a speed of 90 m/s (a moderate-speed movement). The second to fourth scenarios are circular motions at speeds of 100, 500 and 3500 m/s (high-dynamic motions with severe changes in speeds and accelerations). The fifth scenario is a rectangular motion with a speed of about 3200 m/s (highspeed movement with severe accelerations), and the last scenario involves moving in space at a speed of about 7300 m/s (very high-speed movement). The trajectories of the scenarios are shown in Figures 3-8. The speed and acceleration of each scenario are shown in Figures 9-14.
Different modes of movement are used for testing. In the first scenario, there is a variable acceleration of about 5 m/s 2 in the first 100 s, and then acceleration becomes zero. As a result, In second to fourth scenarios, acceleration and velocity are alternating at high rates. In the fifth scenario, due to the rectangular motion, the acceleration at the moments of changing the direction varies suddenly and intensely. In the sixth scenario, the receiver is moving at a very high speed.
The positioning has been performed based on pseudo-range data using LS, RWLS, conventional EKF and the proposed method. The corresponding results are presented below. In Figures 15-20, 3-D root-mean-square errors (RMSEs) of positioning resulting from four previously mentioned methods are shown.
It is observed that in all cases, the LS method has the highest error value. In the RWLS method, the error is reduced but still is high. In addition, the error fluctuations in these two methods are more than the others. The conventional EKF method is more accurate than the two previously mentioned methods; however, in the sudden acceleration case, this method also suffers from fluctuations and instability of the error curve (as in Figure 19). It can be seen that in all cases, the proposed method has the lowest error value and consequently the highest accuracy. It also has a smoother error curve than other methods, and the error is almost constant from the beginning to the end of the simulation. Besides, when the error curve of other methods jumps due to sudden changes in speed or acceleration, the proposed method remains constant and does not increase over time.
The RMSE values of the positioning for different methods in (m) are exhibited in Table 1.
It can be seen that the proposed method has the least positioning error among the methods and has improved the error by more than 70% when compared to the RWLS and conventional EKF methods. Another advantage of the proposed method is that unlike the LS method, it does not depend  Table 2 shows CPU times of methods for each scenario in milliseconds. It can be seen that the CPU time of the proposed method not only is not increased but also is improved more than other methods except the LS method. As mentioned earlier, the advantage of the LS method is its simplicity and therefore has less runtime than other methods. As an advantage, the proposed method does To compare the computational cost of the proposed method and the conventional EKF algorithm, the calculation amounts of EKF and SEKF are presented in Table 3. Details of these calculations are given in the appendices section. In Table 3, l is the observation dimension and n is the state dimension. It is observed that the conventional EKF algorithm has a higher computational cost than the SEKF algorithm due to the two reasons mentioned in Section 6. To provide more clarity on the issue, consider the conditions that existed in the 'Rectangular' data. In this case, n = 11, l = 10 and therefore computational complexity of the conventional EKF algorithm and the SEKF algorithm will be about 10,340 and 5290 flops, respectively.
To further evaluate the performance of the proposed method, the cumulative distribution function (CDF) and probability density function (PDF) of positioning errors are used. The CDF and PDF of four different methods based on
It is shown in Figure 22 that the PDF of the positioning error of the proposed method is more concentrated. Especially most position errors appear within 92 cm and only a few position errors are above 1 m. The error points of other methods are more dispersed. In addition, the position errors of the proposed method appear within 4 m while EKF, RLS, and LS appear within 8, 11, and 28 m, respectively. From the results of CDF and PDF, it can be seen that the proposed method has a better performance than other positioning methods with higher accuracy and more concentrated position errors.
For 'Rectangular' data, some statistical indicators are shown in Table 4. The proposed method exceeds other methods on all indicators, especially on the average accuracy and error (95%) indicators, which goes far beyond others.

| CONCLUSION
The ever-expanding use of GPS constitutes a demand to improve its performance. As the receiver speed increases, its positioning accuracy decreases. The main methods of GPS positioning are the LS method and the KF method. Despite numerous efforts to improve these methods, their accuracy at high speeds or high dynamics is not very good. Here, a comprehensive method is proposed to increase the accuracy of high-speed single-frequency GPS receiver positioning by more than 70% without increasing the computational complexity. Appropriate receiver movement modelling, using the SEKF, correct determination of the Q matrix, and selecting the proper variance estimation method are the basics of the proposed approach. Besides, the proposed method is better than other methods in terms of CPU time. Another advantage of this method is that it does not depend on the starting point. Simulations on different scenarios show that the proposed method performs well in a wide range of speed and acceleration of moving objects. As a suggestion for future study, further models for high-dynamic movement can be explored. The use of UKF as an estimator is another issue that can be examined.