Research on the Error of Global Positioning System Based on Time Series Analysis

Due to the poor dynamic positioning precision of the Global Positioning System (GPS), Time Series Analysis (TSA) and Kalman filter technology are used to construct the positioning error of GPS. According to the statistical characteristics of the autocorrelation function and partial autocorrelation function of sample data, the Autoregressive (AR) model which is based on a Kalman filter is determined, and the error model of GPS is combined with a Kalman filter to eliminate the random error in GPS dynamic positioning data. The least square method is used for model parameter estimation and adaptability tests, and the experimental results show that the absolute value of the maximum error of longitude and latitude, the mean square error of longitude and latitude and average absolute error of longitude and latitude are all reduced, and the dynamic positioning precision after correction has been significantly improved.


Introductions
The United States was the first country to start the new generation of space satellite navigation and position system. Subsequently, Russia and China have studied global navigation positioning systems and developed them rapidly. The research and development of GPS was to serve military needs and, as early as 1972, the United States Navy transit project and the Air Force 621B project realized the application of GPS. With the development of economic globalization, smart phones, vehicles and ships, which are based on the combination of GPS and modern technology, have gradually penetrated into people's daily life [1]. A country's important national infrastructure needs a perfect and stable Global Positioning System, and the precision is an important parameter in the development of GPS [2].
At present, there are two methods to improve the positioning precision of GPS. The first is to use the Differential Global Positioning System (DGPS) [3]. However, this method has several limitations: (1) the equipment limitations: the receiver of a differential signal must be used in practical applications. (2) Restrictions of work area: the work area is limited by the differential network of a wide area. (3) Non-autonomous: the transmitting source is required to improve precision [4]. The other is to use the positioning data which are collected by the receiver for error correction. This method has poor real-time performance and is vulnerable to external interference [5,6]. Therefore, it is particularly important to improve the positioning precision of GPS without changing the hardware device.
Scholars have built the error model of GPS to improve the positioning precision of GPS. In Ref. [7], Zhao Shan et al. used the ARMA (3,2) model to obtain the single error model of a user by the combination of position operation and clock displacement filtering,

The Error Model of GPS
The error data of GPS are a discrete random variable of time series, which is different from the analyzed dynamic data because the time series is a realization of random processes and has a different physical background [14,15]. The processing method of the corresponding model is an approximate description according to the data characteristics, to determine the type of model suitable for the time series [16].
There are three important models of finite order linear in TSA, Moving Average (MA) model, Autoregressive (AR) model, and Autoregressive Moving Average (ARMA) model [17]. The type of model can be determined by analyzing the autocorrelation and partial correlation of the error signal, and the identification methods of the three models are shown in Table 1. As shown in Table 1, the autocorrelation function of MA is truncated, while the partial autocorrelation function is trailing. The autocorrelation function and partial autocorrelation function of ARMA are trailing. The autocorrelation function of AR has the trailing property, while the partial autocorrelation function has the truncation property, where truncation means that the time-order autocorrelation function or partial autocorrelation function is 0 when the order is greater than a constant K. The trailing property means that the autocorrelation function or partial autocorrelation function fluctuates near zero after a certain order.
If a time series is generated by a certain type of model, it should theoretically have corresponding statistical characteristics [18,19]. Therefore, the sample autocorrelation function and partial autocorrelation function of time series can be calculated, and the characteristics can be compared with the characteristics of a theoretical autocorrelation function and partial autocorrelation function of different types of series, and then the model type suitable for the series can be judged [20]. The estimated valueρ k of the autocorrelation function is a measurement to describe the dependence between values of random processes at different times.
After obtaining the error observation data X i (i = 1,2,· · · ,N), the sample autocorrelation functionρ k of the error sequence isρ k =γ k /γ 0 (1) where the estimate value of autocovariance isγ k = 1 N The estimate data of mathematical expectation are X M . Using a numerical method and MATLAB programming, the response simulations of the autocorrelation function and partial autocorrelation function of longitude and latitude are obtained.
By usingγ k instead of γ k , the estimation of the partial autocorrelation function can be obtained recursively. If the autocorrelation functionφ kk of the sample is truncated in step, it can be determined as an AR (M) sequence. Ifφ kk is not truncated, it is an ARMA sequence. From Figures 1a and 2a, the error autocorrelation sequences of longitude and latitude have coordinates fluctuations near zero after order 4, the red line represents the autocorrelation function data of latitude and longitude error. From Figures 1b and 2b, the autocorrelation function of longitude and latitude shows trailing, and the partial autocorrelation function shows truncation, the blue line represents the partial correlation function error data of latitude and longitude error, * represents the data point of latitude and longitude error per second. Therefore, the error of longitude and latitude can be expressed by the AR (M) model, that is, The recursive formula ofρ k can be used to findφ kk by usingρ k instead ofγ k .
tion function k ρ of the error sequence is k k 0ρ γ γ = (1) where the estimate value of autocovariance is The estimate data of mathematical expectation are M X . Using a numerical method and MATLAB programming, the response simulations of the autocorrelation function and partial autocorrelation function of longitude and latitude are obtained.
By using k γ instead of k γ , the estimation of the partial autocorrelation function can be obtained recursively. If the autocorrelation function kk ϕ of the sample is truncated in step, it can be determined as an AR (M) sequence. If kk ϕ is not truncated, it is an ARMA sequence. From Figures 1a and 2a, the error autocorrelation sequences of longitude and latitude have coordinates fluctuations near zero after order 4, the red line represents the autocorrelation function data of latitude and longitude error. From Figures 1b and 2b, the autocorrelation function of longitude and latitude shows trailing, and the partial autocorrelation function shows truncation, the blue line represents the partial correlation function error data of latitude and longitude error, * represents the data point of latitude and longitude error per second. Therefore, the error of longitude and latitude can be expressed by the AR (M) model, that is, The recursive formula of k ρ can be used to find kk ϕ by using k ρ instead of k γ .
Sample autocorrelation of longitude

Preamble Data Processing
We use Google Earth, an interactive electronic map, to capture longitude and latitude information and motion track. From Figure 3a,b, Figure 3a is the Google Earth motion map, a red cross is the trajectory point and a square is an ordinary landmark. Figure 3b is partial latitude and longitude data.  The longitude and latitude data are established in the WGS-84 coordinate system, and the unit is degree-minute second. However, the data in the unit are required, and coordinate projection is required [21]. The WGS-84 coordinate value is projected into the rectangular coordinate value of the Gaussian projection plane by ARCGIS software [22]. As shown in Table 2, the X coordinate is north latitude and the Y coordinate is east longitude.

Preamble Data Processing
We use Google Earth, an interactive electronic map, to capture longitude and latitude information and motion track. From Figure 3a,b, Figure 3a is the Google Earth motion map, a red cross is the trajectory point and a square is an ordinary landmark. Figure 3b is partial latitude and longitude data.

Preamble Data Processing
We use Google Earth, an interactive electronic map, to capture longitude and latitude information and motion track. From Figure 3a,b, Figure 3a is the Google Earth motion map, a red cross is the trajectory point and a square is an ordinary landmark. Figure 3b is partial latitude and longitude data. The longitude and latitude data are established in the WGS-84 coordinate system, and the unit is degree-minute second. However, the data in the unit are required, and coordinate projection is required [21]. The WGS-84 coordinate value is projected into the rectangular coordinate value of the Gaussian projection plane by ARCGIS software [22]. As shown in Table 2, the X coordinate is north latitude and the Y coordinate is east longitude.  The longitude and latitude data are established in the WGS-84 coordinate system, and the unit is degree-minute second. However, the data in the unit are required, and coordinate projection is required [21]. The WGS-84 coordinate value is projected into the rectangular coordinate value of the Gaussian projection plane by ARCGIS software [22]. As shown in Table 2, the X coordinate is north latitude and the Y coordinate is east longitude.
Due to the randomness of the point position in the dynamic environment, the truth value of the observed quantity cannot be obtained directly, which makes the error separation relatively difficult [23]. In order to obtain the position data in the dynamic environment, we adopt the polynomial fitting method to separate the error and select the power series fitting, and the fitting value is the true value of the motion trajectory [24]. According to the polynomial approximation theory, the fitting value compares with the position data continuously collected on this motion trajectory to obtain the estimation of the GPS measurement error value, so as to separate the error data of GPS in the dynamic environment. The estimated value of GPS measurement error is as follows: b i t n+1−i is an n-order power function fitting polynomial model. Coefficient b i is a χ 2 solution with minimum quantity, , ∆x i is the deviation between the original data and the fitted value. If the order of the model polynomial is too low, the fitting is rough. If the order is too high, overfitting will make the model contain data noise. In practical applications, it is necessary to judge whether the fitting is appropriate by the value of fitting factor Q. If Q is close to 0.5, it indicates that the fitting is good.
where P is the probability, N is the number of samples, and n is the order of the model. The calculated Q value is 0.408 after fitting. In Figure 4, there is the large fluctuation of some points which is caused by large-scale movement, not a bad value. In Figure 5, it is the result of error separation of the dynamic position data. Taking the fitting value as the true value estimation of the motion trajectory, and comparing it with the data of position, the error value of dynamic measurements can be obtained by the motion trajectory.

Determination and Parameter Estimation of AR Model Order
It is known that the order of the AR (M) model is the determination of the AR model order when the probability of the partial autocorrelation function falling within the interval 2 √ n exceeds 95% [25]. We take 2 and latitude are truncated in 10 steps. The error model of GPS can be expressed as AR (10), that is, where ϕ i is an autoregressive coefficient, b(k) is the white noise with mean value of 0 and variance of θ 2 .

Determination and Parameter Estimation of AR Model Order
It is known that the order of the AR (M) model is the determination of the AR model order when the probability of the partial autocorrelation function falling within the interval 2 n exceeds 95% [25]. We take

Determination and Parameter Estimation of AR Model Order
It is known that the order of the AR (M) model is the determination of the AR model order when the probability of the partial autocorrelation function falling within the interval 2 n exceeds 95% [25]. We take  Table 3. Partial correlation function of longitude and latitude. As the least square estimation is a precise estimation of model parameters, the estimation precision is high [26]. Therefore, the least square method is used to estimate parameters. It isφ The variance of b(k) is: The autoregressive coefficient and noise variance obtained by the least square method are shown in Table 4.

Applicability Test of AR Model
The applicability of the model is essentially to test whether it is a white noise sequence, and the most important problem is to test the independence of the sequence [27]. When the number of samples is enough, the autocorrelation functions of residual error are uncorrelated and approximate to the normal distribution.
The green symbol is the residual value within the confidence interval, and the red symbol is the residual value outside the confidence interval. The steadiness of the model is also an indicator to test the residual sequence. It can be seen from Figure 6 that the autocorrelation coefficient of residual sequence fluctuates randomly with 0 as the mean in the 95% confidence boundary, and the residual sequence data of the longitude and latitude are stable, which indicates that the model has a high degree of fitting with the actual system and the model meets the requirements. Therefore, the AR (10) model is more suitable for GPS longitude and latitude error time series. The green symbol is the residual value within the confidence interval, and the red symbol is the residual value outside the confidence interval. The steadiness of the model is also an indicator to test the residual sequence. It can be seen from Figure 6 that the autocorrelation coefficient of residual sequence fluctuates randomly with 0 as the mean in the 95% confidence boundary, and the residual sequence data of the longitude and latitude are stable, which indicates that the model has a high degree of fitting with the actual system and the model meets the requirements. Therefore, the AR (10) model is more suitable for GPS longitude and latitude error time series.

Discretization of State Equations of Continuous Systems
The actual physical system is generally continuous, and the dynamic characteristics are described by continuous differential equations. Therefore, the discretizations of the system equation and the observation equation are needed [28,29].
The system state equation describing the dynamic characteristics of the physical system is: where the driving source v(t) of the system is the white noise process, which is where q is v(t) variance intensity matrix. σ(t-τ) is a function of Dirac σ.
According to the linear system theory, the discretization of the system state equation is: where the one-step transfer matrix a k+1 k (t ,t ) satisfies the equation: where k k k+1 k f =f(t ),T=t -t .The equation is the real-time calculation formula of a one-step transfer matrix.
The discretization state equation of continuous system also includes the equivalent discretization of the excited white noise process v(t) .

Discretization of State Equations of Continuous Systems
The actual physical system is generally continuous, and the dynamic characteristics are described by continuous differential equations. Therefore, the discretizations of the system equation and the observation equation are needed [28,29].
The system state equation describing the dynamic characteristics of the physical system is: .
where the driving source v(t) of the system is the white noise process, which is where q is v(t) variance intensity matrix. σ(t-τ) is a function of Dirac σ. According to the linear system theory, the discretization of the system state equation is: where the one-step transfer matrix a (t k+1 , t k satisfies the equation: where f k = f(t k ), T = t k+1 −t k . The equation is the real-time calculation formula of a onestep transfer matrix.
The discretization state equation of continuous system also includes the equivalent discretization of the excited white noise process v(t). Equation (11) can be abbreviated as: (14) where X k = X(t k ). Then, for V(k) defined in Equation (13), it is: (15) where, δ kj is Kronecker δ function. The variance matrix Q k of V k satisfies the following equation: where, g k = g(t k ), the Equation (16) is the real-time calculation formula of Q k .

The Basic Equation of Discrete Kalman Filter
Kalman is linear minimum variance estimation [30]. For Kalman model, the state equation and observation equation of discrete linear system are respectively: where, a k,k−1 is one step transition matrix of time, which is t k−1 to t k . c k is the measurement matrix. ξ k−1 is the system noise driving. W k is the observation noise sequence. V k−1 is the excitation noise sequence of the system. Both V k and W k are satisfied where, Q k is the variance matrix of the system noise sequence, it is a non-negative matrix. R k is the variance matrix of noise sequence on both sides, it is positive definite matrix. δ kj is Kronecker δ function. If X k is conformed to Equation (19), the measured value Y k is conformed to Equation (20), the system noise V k and the measurement noise W k are conformed to Equation (21), and the system noise variance matrix Q k is non-negative definite, the measurement noise variance matrix R K is positive definite, and the measurement of time K is Y k , the equation can be solved for X K estimationX K . Then where, X k|k−1 is the one-step prediction equation of state. P k|k−1 is the one-step prediction of mean square error. K k is the gain equation of filter.X k is the estimation equation of filter. p k is the optimal mean square error at time K. As long as the initial valueX 0 and P 0 is given, according to the measurement Y k at time k, the estimated stateX K (k = 1, 2, 3, · · · ) at time k can be deduced.

Kalman Filter Based on AR Model
The state vector of GPS is . . .
The model is based on x(k) = 10 ∑ n=1 ϕ n x(k − n) + b(k), the state space model of Kalman filter model is Equations (9) and (10), where V k and W k have statistical characteristics of Equation (21). The process noise turn the difference equation of AR (10) model into state Equation (24), it is . . .
. . . (24) can also be abbreviated as The observation equation is where, C = 0 0 · · · 0 1 . Because the statistical characteristics of V(k) and W(k) are consistent with Kalman filter, the recursive expression of Kalman filter based on AR model is show as: where, X k|k−1 is the one-step prediction equation of filter state. P k|k−1 is the one-step prediction of mean square error. K k is the gain equation of filter. If the Kalman filter gain value K k is very small, the filtering result is closer to the recursive result of the system state estimation value. If K k is large, the filtering result is closer to the state variable calculated of the observed value.X k is the estimation equation of filter. p k is the optimal mean square error at time K. There are two ways to deal with Q: one is that Q is a certain value. The second is that Q is an uncertain random variable. Therefore, the Q value in this paper is a definite value, which is the process noise variance of the system, and its value is σ * I 10order . When the state transition process has been determined, the smaller of Q is the better. When Q gradually increases, the convergence of filter slows down and the disturbance of the state variable becomes larger. The value of R is related to the characteristics of the device and is the input value of the filter.

AR (10) Model Kalman Filter Experiment
The dynamic error model (AR model) of GPS uses a Kalman filter, Q/R is 0.1/0.5, and the order of the AR model is determined to be 10. The error filtering results of longitude and latitude are shown in Figure 7.
culated of the observed value. k X is the estimation equation of filter. k p is the optimal mean square error at time K.
There are two ways to deal with Q: one is that Q is a certain value. The second is that Q is an uncertain random variable. Therefore, the Q value in this paper is a definite value, which is the process noise variance of the system, and its value is 10order I σ * . When the state transition process has been determined, the smaller of Q is the better. When Q gradually increases, the convergence of filter slows down and the disturbance of the state variable becomes larger. The value of R is related to the characteristics of the device and is the input value of the filter.

AR (10) Model Kalman Filter Experiment
The dynamic error model (AR model) of GPS uses a Kalman filter, Q/R is 0.1/0.5, and the order of the AR model is determined to be 10. The error filtering results of longitude and latitude are shown in Figure 7. In Figure 7, the blue curve is the error data of original longitude and latitude, the green curve is the error of longitude and latitude after the Kalman filter, and the red curve is the ideal error. The blue and green lines do not completely coincide, and the green lines after the Kalman filter are smoother. There are several points with large amplitude in the figure, which are caused by excessive movement when collecting data. In order to analyze Figure 7, Table 5 is shown below for explanation. In Figure 7, the blue curve is the error data of original longitude and latitude, the green curve is the error of longitude and latitude after the Kalman filter, and the red curve is the ideal error. The blue and green lines do not completely coincide, and the green lines after the Kalman filter are smoother. There are several points with large amplitude in the figure, which are caused by excessive movement when collecting data. In order to analyze Figure 7, Table 5 is shown below for explanation. The analysis of data is carried out, and the absolute value of maximum error, mean square error, and mean absolute error are selected to reflect the changes in data before and after the filtering model. It can be seen from Table 5 that the absolute value of the maximum error of longitude and latitude after filtering is reduced by 30.50% and 37.26%. The mean square errors of longitude and latitude were reduced by 8.878% and 10.12%. The mean absolute errors of longitude and latitude were also reduced by 5.695% and 5.336%.

Compared Experiment of AR Model
In order to further prove the effectiveness of the model in this paper, a new set of GPS data is collected, and an AR model based on a different order Kalman filter is established and compared with the model. Q/R is 0.1/0.5, the longitude and latitude error curve of each model is shown in Figure 8.

Compared Experiment of AR Model
In order to further prove the effectiveness of the model in this paper, a new set of GPS data is collected, and an AR model based on a different order Kalman filter is established and compared with the model. Q/R is 0.1/0.5, the longitude and latitude error curve of each model is shown in Figure 8. In Figure 8, the red line is the error of expectation, the blue line is the error of measurement, and the brown line the error of the AR (10) model. In order to more intuitively compare the error of each model, the longitude and latitude error data of each model are given in Tables 6 and 7.  In Figure 8, the red line is the error of expectation, the blue line is the error of measurement, and the brown line the error of the AR (10) model. In order to more intuitively compare the error of each model, the longitude and latitude error data of each model are given in Tables 6 and 7. It can be seen from Tables 6 and 7 that the absolute value of the maximum error of longitude and latitude in this paper has been reduced by 41.67% and 52.25%, respectively. The mean square errors of longitude and latitude were reduced by 37.41% and 39.48%, respectively. The mean absolute errors of longitude and latitude were also reduced by 36.53% and 23.26%, respectively. Although the errors of other order models are also greatly reduced with a Kalman filter, the AR (10) model in this paper is selected after order determination and adaptability tests, so the positioning accuracy of this model is better than other models.

Conclusions
In this study, the error of GPS is researched, and the statistical characteristics of GPS are analyzed and simulated. The least square method is applied to estimate the parameters, and obtain the basic equation of a discrete Kalman filter by the continuous Kalman filter. To eliminate the random error of GPS dynamic positioning data, the error model of GPS is combined with a Kalman filter and the experimental results show that the smaller the mean square error, the better the precision of the prediction model on the experimental data. The dynamic positioning precision after correction has been significantly improved, and this method can effectively improve the positioning precision of GPS and support for the application of GPS.