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

Abstract

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.

1. 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, which has a certain stability. In Ref. [8], Liu Di et al. made a long-term observation of a static point to obtain the same characteristics of error sequence and elevation error sequence between longitude and latitude, and finally established the AR (n) error model of GPS. In Ref. [9], Wang Rong et al. analyzed the error of GPS and established the error model of GPS. In Ref. [10], Zhiqiang Liu et al. modeled the time-selective channel as an AR process and used a Kalman filter to track the time change. In Ref. [11], Tan-Jan Ho et al. proposed a framework modeling which is based on multiple AR models, and developed the channel predictor MAR. In Ref. [12], Christos Komnanakis et al. used a Kalman filter to track a low-order autoregressive model, which is similar to the change in multiple input and multiple output channels.

The models used AR models based on time series by the scholars. If the output sequence of a model is ${\mathrm{x}}_{\mathrm{t}}$, there is a non-linear problem of parameter estimation, and the algorithm is complex and difficult. In order to solve the non-linear problems and effectively improve the positioning precision of GPS, we use the TAS and Kalman filter to analyze the characteristics of positioning error data, set up the Autoregressive (AR) model of GPS, complete the parameter estimation of the AR model by the least square method, and estimate the state of the system by a Kalman filter, so as to make the positioning precision more reliable and accurate [13].

2. 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 ${\widehat{\rho }}_{\mathrm{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 ${\mathrm{X}}_{\mathrm{i}}\mathrm{\left(}\mathrm{i}\mathrm{=}\mathrm{1}\mathrm{,}\mathrm{2}\mathrm{,}\mathrm{\cdots }\mathrm{,}\mathrm{N}\mathrm{\right)}$, the sample autocorrelation function ${\widehat{\rho }}_{\mathrm{k}}$ of the error sequence is

(1)${\widehat{\rho }}_{\mathrm{k}}={\widehat{\gamma }}_{\mathrm{k}}/{\mathrm{\gamma }}_0$

where the estimate value of autocovariance is ${\widehat{\mathrm{\gamma }}}_{\mathrm{k}}=\frac{1}{\mathrm{N}}{\displaystyle \sum _{\mathrm{i}=\mathrm{1}}^{\mathrm{N}-\mathrm{k}}\left({\mathrm{X}}_{\mathrm{i}}-{\mathrm{X}}_{\mathrm{M}}\right)\left({\mathrm{X}}_{\mathrm{i}+\mathrm{k}}-{\mathrm{X}}_{\mathrm{M}}\right)}$, $\mathrm{k}=0,1,\cdots$. The estimate data of mathematical expectation are ${\mathrm{X}}_{\mathrm{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 ${\widehat{\mathrm{\gamma }}}_{\mathrm{k}}$ instead of ${\mathrm{\gamma }}_{\mathrm{k}}$, the estimation of the partial autocorrelation function can be obtained recursively. If the autocorrelation function ${\widehat{\phi }}_{\mathrm{k}\mathrm{k}}$ of the sample is truncated in step, it can be determined as an AR (M) sequence. If ${\widehat{\phi }}_{\mathrm{k}\mathrm{k}}$ is not truncated, it is an ARMA sequence. From Figure 1a and Figure 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 Figure 1b and Figure 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,

(2)$\{\begin{array}{c}{\widehat{\phi }}_{11}={\widehat{\mathrm{\gamma }}}_1/{\widehat{\mathrm{\gamma }}}_0\\ {\widehat{\phi }}_{\begin{array}{cc}\mathrm{k}+1& \mathrm{k}+1\end{array}}=\left({\widehat{\mathrm{\gamma }}}_{\mathrm{k}+1}-{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{k}}{\widehat{\phi }}_{\mathrm{k}\mathrm{j}}{\widehat{\mathrm{\gamma }}}_{\mathrm{k}+1-\mathrm{j}}}\right)/\left({\widehat{\mathrm{\gamma }}}_0-{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{k}}{\widehat{\phi }}_{\mathrm{k}\mathrm{j}}{\widehat{\mathrm{\gamma }}}_{\mathrm{j}}}\right)\\ {\widehat{\phi }}_{\begin{array}{cc}\mathrm{k}+1& \mathrm{j}\end{array}}={\widehat{\phi }}_{\mathrm{k}\mathrm{j}}-{\widehat{\phi }}_{\begin{array}{cc}\mathrm{k}+1& \mathrm{k}+1\end{array}}{\widehat{\phi }}_{\begin{array}{cc}\mathrm{k}& \mathrm{k}+1-\mathrm{j}\end{array}}\mathrm{,}\mathrm{j}\mathrm{=}\mathrm{1}\mathrm{,}\mathrm{2}\mathrm{,}\mathrm{\cdots }\mathrm{,}\mathrm{k}\end{array}.$

The recursive formula of ${\widehat{\rho }}_{\mathrm{k}}$ can be used to find ${\widehat{\phi }}_{\mathrm{k}\mathrm{k}}$ by using ${\widehat{\rho }}_{\mathrm{k}}$ instead of ${\widehat{\mathrm{\gamma }}}_{\mathrm{k}}$.

3. The Establishment of the Error Model of GPS

3.1. 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.

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:

(3)${\widehat{\mathrm{X}}}_{\mathrm{G}\mathrm{P}\mathrm{S}}\mathrm{=}{\mathrm{x}}_{\mathrm{i}}\mathrm{-}\mathrm{f}\mathrm{\left(}\mathrm{b}\mathrm{,}\mathrm{t}\mathrm{\right)}$

where ${\mathrm{x}}_{\mathrm{i}}$ is the measured value. $\mathrm{f}\left(\mathrm{b},\mathrm{t}\right)={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}+1}{\mathrm{b}}_{\mathrm{i}}{\mathrm{t}}^{\mathrm{n}+1-\mathrm{i}}}$ is an n-order power function fitting polynomial model. Coefficient ${\mathrm{b}}_{\mathrm{i}}$ is a ${\mathrm{\chi }}^2$ solution with minimum quantity, ${\mathrm{\chi }}^2\left(\mathrm{b}\right)={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{N}}{\left(\frac{{\mathrm{x}}_{\mathrm{i}}-\mathrm{f}\left(\mathrm{b},{\mathrm{t}}_{\mathrm{i}}\right)}{{\mathsf{\Delta }\mathrm{x}}_{\mathrm{i}}}\right)}^2$, ${\mathsf{\Delta }\mathrm{x}}_{\mathrm{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.

(4)$\mathrm{Q}\left({\mathrm{\chi }}^2,\mathrm{N}-\mathrm{n}-1\right)=1-\mathrm{P}({\mathrm{\chi }}^2<\left(\mathrm{N}-\mathrm{n}-1\right))$

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.

3.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 $\left|\frac{2}{\sqrt{\mathrm{n}}}\right|$ exceeds 95% [25]. We take $\frac{2}{\sqrt{\mathrm{n}}}=\frac{2}{\sqrt{329}}=0.11$ and obtain Table 3 and Table 4 by a MATLAB programming calculation, and the partial autocorrelation functions of longitude and latitude are truncated in 10 steps. The error model of GPS can be expressed as AR (10), that is,

(5)$\mathrm{y}\left(\mathrm{k}\right)={\displaystyle \sum _{\mathrm{i}=1}^{10}{\phi }_{\mathrm{i}}\mathrm{y}\left(\mathrm{k}-\mathrm{i}\right)+\mathrm{b}\left(\mathrm{k}\right)}$

where ${\phi }_{\mathrm{i}}$ is an autoregressive coefficient, $\mathrm{b}\left(\mathrm{k}\right)$ is the white noise with mean value of 0 and variance of ${\mathrm{\theta }}^2$.

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

(6)${\widehat{\phi }}_{\mathrm{L}}={\left({\mathrm{X}}^{\mathrm{T}}\mathrm{X}\right)}^{-1}{\mathrm{X}}^{\mathrm{T}}\mathrm{Y}.$

The variance of $\mathrm{b}(\mathrm{k})$ is:

(7)${\mathrm{\theta }}^2=\frac{1}{\mathrm{I}-\mathrm{n}}{\displaystyle \sum _{\mathrm{t}=\mathrm{n}+1}^{\mathrm{I}}\le \left(\mathrm{y}\left(\mathrm{k}\right)-{\phi }_1\mathrm{y}\left(\mathrm{k}-1\right)-\cdots -{\phi }_{\mathrm{i}}\mathrm{y}\left({\mathrm{k}}^2-\mathrm{n}\right)\right)}.$

The autoregressive coefficient and noise variance obtained by the least square method are shown in Table 4.

3.3. 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.

(8)${\widehat{\rho }}_{\mathrm{k}}={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{i}-\mathrm{k}}{\widehat{\mathrm{b}}}_1{\widehat{\mathrm{b}}}_{1+\mathrm{k}}/\left[{\displaystyle \sum _{\mathrm{t}-1}^{\mathrm{i}}{\widehat{\mathrm{b}}}^2}\right]}.$

At the significance level $\alpha =0.05$, $\left|{\widehat{\rho }}_{\mathrm{k}}\right|\le \frac{1.96}{{\mathrm{N}}^{0.5}}$, it is acceptable to assume that ${\widehat{\rho }}_{\mathrm{k}}=0$ is independent.

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.

4. The Application of Kalman Filter in the Error Model

4.1. 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:

(9)$\dot{\mathrm{X}}(\mathrm{t})=\mathrm{f}(\mathrm{t})\mathrm{X}(\mathrm{t})+\mathrm{g}(\mathrm{t})\mathrm{v}(\mathrm{t})$

where the driving source $\mathrm{v}\left(\mathrm{t}\right)$ of the system is the white noise process, which is

(10)$\{\begin{array}{c}\mathrm{E}[\mathrm{v}(\mathrm{t})]=0\\ \mathrm{E}[\mathrm{v}(\mathrm{t}){\mathrm{v}}^{\mathrm{T}}(\tau )]=\mathrm{q}\sigma (\mathrm{t}-\tau )\end{array}$

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:

(11)$\mathrm{X}({\mathrm{t}}_{\mathrm{k}+1})=\mathrm{a}({\mathrm{t}}_{\mathrm{k}+1},{\mathrm{t}}_{\mathrm{k}})\mathrm{X}({\mathrm{t}}_{\mathrm{k}})+{\displaystyle {\mathrm{\int }}_{{\mathrm{t}}_{\mathrm{k}}}^{{\mathrm{t}}_{\mathrm{k}+1}}\mathrm{a}({\mathrm{t}}_{\mathrm{k}+1},\tau )}\mathrm{g}(\tau )\mathrm{v}(\tau )\mathrm{d}\tau$

where the one-step transfer matrix a ${(\mathrm{t}}_{\mathrm{k}+1}{,\mathrm{t}}_{\mathrm{k}})$ satisfies the equation:

(12)${\mathrm{a}}_{\mathrm{k}+1,\mathrm{k}}=\mathrm{a}({\mathrm{t}}_{\mathrm{k}+1},{\mathrm{t}}_{\mathrm{k}})=\mathrm{I}+\mathrm{T}{\mathrm{f}}_{\mathrm{k}}+\frac{{\mathrm{T}}^2}{2!}{\mathrm{f}}_{\mathrm{k}}^2+\frac{{\mathrm{T}}^3}{3!}{\mathrm{f}}_{\mathrm{k}}^3+\frac{{\mathrm{T}}^4}{4!}{\mathrm{f}}_{\mathrm{k}}^4+\cdots$

where ${\mathrm{f}}_{\mathrm{k}}{=\mathrm{f}(\mathrm{t}}_{\mathrm{k}}{),\mathrm{T}=\mathrm{t}}_{\mathrm{k}+1}{-\mathrm{t}}_{\mathrm{k}}$. 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 $\mathrm{v}\left(\mathrm{t}\right)$.

(13)${\mathrm{V}}_{\mathrm{k}}={\displaystyle {\int }_{{\mathrm{t}}_{\mathrm{k}}}^{{\mathrm{t}}_{\mathrm{k}+1}}\mathrm{\varphi }({\mathrm{t}}_{\mathrm{k}+1},\tau )\mathrm{g}(\tau )\mathrm{v}(\tau )\mathrm{d}\tau }.$

Equation (11) can be abbreviated as:

(14)${\mathrm{X}}_{\mathrm{k}+1}={\mathrm{a}}_{\mathrm{k}+1,\mathrm{k}}{\mathrm{X}}_{\mathrm{k}}+{\mathrm{V}}_{\mathrm{k}}$

where ${\mathrm{X}}_{\mathrm{k}}{=\mathrm{X}(\mathrm{t}}_{\mathrm{k}})$. Then, for $\mathrm{V}\left(\mathrm{k}\right)$ defined in Equation (13), it is:

(15)$\{\begin{array}{c}\mathrm{E}[{\mathrm{V}}_{\mathrm{k}}]=0\\ \mathrm{E}[{\mathrm{V}}_{\mathrm{k}}{\mathrm{V}}_{\mathrm{j}}^{\mathrm{T}}]={\displaystyle {\int }_{{\mathrm{t}}_{\mathrm{k}}}^{{\mathrm{t}}_{\mathrm{k}+1}}\mathrm{\varphi }({\mathrm{t}}_{\mathrm{k}+1},{\mathrm{t}}_{\mathrm{k}})\mathrm{g}(\mathrm{t})\mathrm{q}{\mathrm{g}}^{\mathrm{T}}(\mathrm{t}){\mathrm{\varphi }}^{\mathrm{T}}({\mathrm{t}}_{\mathrm{k}+1},{\mathrm{t}}_{\mathrm{k}})\mathrm{d}\mathrm{t}=\mathrm{Q}{\mathsf{\delta }}_{\mathrm{k}\mathrm{j}}}\end{array}$

where, ${\mathsf{\delta }}_{\mathrm{kj}}$ is Kronecker $\mathsf{\delta }$ function. The variance matrix ${\mathrm{Q}}_{\mathrm{k}}$ of ${\mathrm{V}}_{\mathrm{k}}$ satisfies the following equation:

(16)${\mathrm{M}}_1={\mathrm{g}}_{\mathrm{k}}\mathrm{q}{\mathrm{g}}_{\mathrm{k}}{}^{\mathrm{T}}$

(17)${\mathrm{M}}_{\mathrm{i}+1}={\mathrm{f}}_{\mathrm{k}}\mathrm{q}{\mathrm{g}}_{\mathrm{k}}{}^{\mathrm{T}}$

(18)${\mathrm{Q}}_{\mathrm{k}}=\mathrm{T}{\mathrm{M}}_1+\frac{{\mathrm{T}}^2}{2!}{\mathrm{M}}_2+\frac{{\mathrm{T}}^3}{3!}{\mathrm{M}}_3+\frac{{\mathrm{T}}^4}{4!}{\mathrm{M}}_4+\cdots$

where, ${\mathrm{g}}_{\mathrm{k}}{=\mathrm{g}(\mathrm{t}}_{\mathrm{k}})$, the Equation (16) is the real-time calculation formula of ${\mathrm{Q}}_{\mathrm{k}}$.

4.2. 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:

(19)${\mathrm{X}}_{\mathrm{k}}={\mathrm{a}}_{\mathrm{k},\mathrm{k}-1}{\mathrm{X}}_{\mathrm{k}-1}+{\mathrm{\xi }}_{\mathrm{k}-1}{\mathrm{V}}_{\mathrm{k}-1}$

(20)${\mathrm{Y}}_{\mathrm{k}}={\mathrm{c}}_{\mathrm{k}}{\mathrm{X}}_{\mathrm{k}}+{\mathrm{W}}_{\mathrm{k}}$

where, ${\mathrm{a}}_{\mathrm{k},\mathrm{k}-1}$ is one step transition matrix of time, which is ${\mathrm{t}}_{\mathrm{k}-1}$ to ${\mathrm{t}}_{\mathrm{k}}$. ${\mathrm{c}}_{\mathrm{k}}$ is the measurement matrix. ${\mathsf{\xi }}_{\mathrm{k}-1}$ is the system noise driving. ${\mathrm{W}}_{\mathrm{k}}$ is the observation noise sequence. ${\mathrm{V}}_{\mathrm{k}-1}$ is the excitation noise sequence of the system. Both ${\mathrm{V}}_{\mathrm{k}}$ and ${\mathrm{W}}_{\mathrm{k}}$ are satisfied

(21)$\{\begin{array}{c}{\mathrm{E}(\mathrm{V}}_{\mathrm{K}}{)=0,\mathrm{E}(\mathrm{V}}_{\mathrm{K}}{\mathrm{V}}_{\mathrm{j}}^{\mathrm{T}}{)=0,=\mathrm{Q}}_{\mathrm{K}}{\mathsf{\delta }}_{\mathrm{Kj}}\\ {\mathrm{E}(\mathrm{W}}_{\mathrm{K}}{)=0,\mathrm{E}(\mathrm{W}}_{\mathrm{K}}{\mathrm{W}}_{\mathrm{j}}^{\mathrm{T}}{)=0,=\mathrm{R}}_{\mathrm{K}}{\mathsf{\delta }}_{\mathrm{Kj}}\\ {\mathrm{E}(\mathrm{W}}_{\mathrm{K}}{\mathrm{V}}_{\mathrm{j}}^{\mathrm{T}})=0\end{array}$

where, ${\mathrm{Q}}_{\mathrm{k}}$ is the variance matrix of the system noise sequence, it is a non-negative matrix. ${\mathrm{R}}_{\mathrm{k}}$ is the variance matrix of noise sequence on both sides, it is positive definite matrix. ${\mathsf{\delta }}_{\mathrm{kj}}$ is Kronecker δ function.

If ${\mathrm{X}}_{\mathrm{k}}$ is conformed to Equation (19), the measured value ${\mathrm{Y}}_{\mathrm{k}}$ is conformed to Equation (20), the system noise ${\mathrm{V}}_{\mathrm{k}}$ and the measurement noise ${\mathrm{W}}_{\mathrm{k}}$ are conformed to Equation (21), and the system noise variance matrix ${\mathrm{Q}}_{\mathrm{k}}$ is non-negative definite, the measurement noise variance matrix ${\mathrm{R}}_{\mathrm{K}}$ is positive definite, and the measurement of time K is ${\mathrm{Y}}_{\mathrm{k}}$, the equation can be solved for ${\mathrm{X}}_{\mathrm{K}}$ estimation ${\widehat{\mathrm{X}}}_{\mathrm{K}}$.

Then

(22)$\{\begin{array}{c}{\mathrm{X}}_{\mathrm{k}|\mathrm{k}-1}={\mathrm{a}}_{\mathrm{k},\mathrm{k}-1}{\widehat{\mathrm{X}}}_{\mathrm{k}}\\ {\mathrm{P}}_{\mathrm{k}|\mathrm{k}-1}={\mathrm{a}}_{\mathrm{k},\mathrm{k}-1}{\mathrm{P}}_{\mathrm{k}-1}{\mathrm{a}}^{\mathrm{T}}{}_{\mathrm{k},\mathrm{k}-1}+{\mathrm{\xi }}_{\mathrm{k}-1}{\mathrm{Q}}_{\mathrm{k}-1}{\mathrm{\xi }}_{\mathrm{k}-1}^{\mathrm{T}}\\ {\mathrm{K}}_{\mathrm{k}}={\mathrm{P}}_{\mathrm{k}|\mathrm{k}-1}{\mathrm{c}}_{\mathrm{k}}^{\mathrm{T}}{\left({\mathrm{c}}_{\mathrm{k}}{\mathrm{P}}_{\mathrm{k}|\mathrm{k}-1}{\mathrm{c}}_{\mathrm{k}}^{\mathrm{T}}+{\mathrm{R}}_{\mathrm{k}}\right)}^{-1}\\ {\widehat{\mathrm{X}}}_{\mathrm{k}}={\widehat{\mathrm{X}}}_{\mathrm{k}|\mathrm{k}-1}+{\mathrm{K}}_{\mathrm{k}}\left({\mathrm{Y}}_{\mathrm{k}}-{\mathrm{c}}_{\mathrm{k}}{\widehat{\mathrm{X}}}_{\mathrm{k}|\mathrm{k}-1}\right)\\ {\mathrm{p}}_{\mathrm{k}}=\left(\mathrm{I}-{\mathrm{K}}_{\mathrm{k}}{\mathrm{c}}_{\mathrm{k}}\right){\mathrm{P}}_{\mathrm{k}|\mathrm{k}-1}{\left(\mathrm{I}-{\mathrm{K}}_{\mathrm{k}}{\mathrm{c}}_{\mathrm{k}}\right)}^{\mathrm{T}}+{\mathrm{K}}_{\mathrm{k}}{\mathrm{R}}_{\mathrm{k}}{\mathrm{K}}_{\mathrm{k}}{}^{\mathrm{T}}\end{array}$

where, ${\mathrm{X}}_{\mathrm{k}|\mathrm{k}-1}$ is the one-step prediction equation of state. ${\mathrm{P}}_{\mathrm{k}|\mathrm{k}-1}$ is the one-step prediction of mean square error. ${\mathrm{K}}_{\mathrm{k}}$ is the gain equation of filter. ${\widehat{\mathrm{X}}}_{\mathrm{k}}$ is the estimation equation of filter. ${\mathrm{p}}_{\mathrm{k}}$ is the optimal mean square error at time K.

As long as the initial value ${\widehat{\mathrm{X}}}_0$ and ${\mathrm{P}}_0$ is given, according to the measurement ${\mathrm{Y}}_{\mathrm{k}}$ at time k, the estimated state ${\widehat{\mathrm{X}}}_{\mathrm{K}}(\mathrm{k}=1,2,3,\cdots )$ at time k can be deduced.

4.3. Kalman Filter Based on AR Model

The state vector of GPS is

(23)$\mathrm{x}\left(\mathrm{k}-1\right)=\left[\begin{array}{c}\mathrm{x}\left(\mathrm{k}-10\right)\\ \mathrm{x}\left(\mathrm{k}-9\right)\\ ⋮\\ \mathrm{x}\left(\mathrm{k}-1\right)\end{array}\right]$

The model is based on $\mathrm{x}\left(\mathrm{k}\right)={\displaystyle \sum _{\mathrm{n}=1}^{10}{\phi }_{\mathrm{n}}\mathrm{x}\left(\mathrm{k}-\mathrm{n}\right)+\mathrm{b}\left(\mathrm{k}\right)}$, the state space model of Kalman filter model is Equations (9) and (10), where ${\mathrm{V}}_{\mathrm{k}}$ and ${\mathrm{W}}_{\mathrm{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)$\left[\begin{array}{c}\mathrm{x}\left(\mathrm{k}-9\right)\\ \mathrm{x}\left(\mathrm{k}-8\right)\\ ⋮\\ \mathrm{x}\left(\mathrm{k}\right)\end{array}\right]=\left[\begin{array}{ccccc}0& 1& 0& \cdots & 0\\ 0& 0& 1& 0& 0\\ ⋮& ⋮& & & ⋮\\ 0& 0& \cdots & \cdots & 1\\ {\phi }_{10}& {\phi }_9& \cdots & \cdots & {\phi }_1\end{array}\right]\left[\begin{array}{c}\mathrm{x}\left(\mathrm{k}-10\right)\\ \mathrm{x}\left(\mathrm{k}-9\right)\\ ⋮\\ \mathrm{x}\left(\mathrm{k}-1\right)\end{array}\right]+\left[\begin{array}{c}0\\ 0\\ ⋮\\ 0\\ 1\end{array}\right]\mathrm{V}\left(\mathrm{k}\right)$

Let $\left[\begin{array}{ccccc}0& 1& 0& \cdots & 0\\ 0& 0& 1& 0& 0\\ ⋮& ⋮& & & ⋮\\ 0& 0& \cdots & \cdots & 1\\ {\phi }_{10}& {\phi }_9& \cdots & \cdots & {\phi }_1\end{array}\right]=\mathrm{A},\left[\begin{array}{c}0\\ 0\\ ⋮\\ 0\\ 1\end{array}\right]=\mathrm{B}$, Equation (24) can also be abbreviated as

(25)$\mathrm{x}\left(\mathrm{k}\right)=\mathrm{A}\mathrm{x}(\mathrm{k}-1)+\mathrm{B}\mathrm{V}\left(\mathrm{k}\right)$

The observation equation is

(26)$\mathrm{y}\left(\mathrm{k}\right)=\mathrm{x}\left(\mathrm{k}\right)+\mathrm{W}\left(\mathrm{k}\right)=\mathrm{C}\mathrm{x}\left(\mathrm{k}\right)+\mathrm{W}\left(\mathrm{k}\right)$

where, $\mathrm{C}=\left[\begin{array}{ccccc}0& 0& \cdots & 0& 1\end{array}\right]$.

Because the statistical characteristics of $\mathrm{V}\left(\mathrm{k}\right)$ and $\mathrm{W}\left(\mathrm{k}\right)$ are consistent with Kalman filter, the recursive expression of Kalman filter based on AR model is show as:

(27)$\{\begin{array}{c}{\widehat{\mathrm{X}}}_{\mathrm{k}|\mathrm{k}-1}=\mathrm{A}{\widehat{\mathrm{X}}}_{\mathrm{k}-1}\\ {\mathrm{P}}_{\mathrm{k}|\mathrm{k}-1}=\mathrm{A}{\mathrm{P}}_{\mathrm{k}-1}{\mathrm{A}}^{\mathrm{T}}+\mathrm{B}{\mathrm{Q}}_{\mathrm{k}-1}{\mathrm{B}}^{\mathrm{T}}\\ {\mathrm{K}}_{\mathrm{k}}={\mathrm{P}}_{\mathrm{k}|\mathrm{k}-1}{\mathrm{c}}^{\mathrm{T}}{\left(\mathrm{c}{\mathrm{P}}_{\mathrm{k}|\mathrm{k}-1}{\mathrm{c}}^{\mathrm{T}}+\mathrm{R}\right)}^{-1}\\ {\widehat{\mathrm{X}}}_{\mathrm{k}}={\widehat{\mathrm{X}}}_{\mathrm{k}|\mathrm{k}-1}+{\mathrm{K}}_{\mathrm{k}}\left({\mathrm{Y}}_{\mathrm{k}}-{\mathrm{c}}_{\mathrm{k}}{\widehat{\mathrm{X}}}_{\mathrm{k}|\mathrm{k}-1}\right)\\ {\mathrm{P}}_{\mathrm{k}}=\left(\mathrm{I}-{\mathrm{K}}_{\mathrm{k}}{\mathrm{c}}_{\mathrm{k}}\right){\mathrm{P}}_{\mathrm{k}|\mathrm{k}-1}\end{array}$

where, ${\mathrm{X}}_{\mathrm{k}|\mathrm{k}-1}$ is the one-step prediction equation of filter state. ${\mathrm{P}}_{\mathrm{k}|\mathrm{k}-1}$ is the one-step prediction of mean square error. ${\mathrm{K}}_{\mathrm{k}}$ is the gain equation of filter. If the Kalman filter gain value ${\mathrm{K}}_{\mathrm{k}}$ is very small, the filtering result is closer to the recursive result of the system state estimation value. If ${\mathrm{K}}_{\mathrm{k}}$ is large, the filtering result is closer to the state variable calculated of the observed value. ${\widehat{\mathrm{X}}}_{\mathrm{k}}$ is the estimation equation of filter. ${\mathrm{p}}_{\mathrm{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 $\mathrm{\sigma }\ast {\mathrm{I}}_{10\mathrm{order}}$. 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.

5. Simulation and Analysis

5.1. 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.

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%.

5.2. 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 Table 6 and Table 7.

It can be seen from Table 6 and Table 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.

6. 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.

Author Contributions

Conceptualization, L.S.; Data curation, L.Z., P.X. and Z.L.; Formal analysis, W.Z. and S.L. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The raw/processed data required to reproduce these findings cannot be shared at this time as the data also form part of an ongoing study.

Conflicts of Interest

We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Autocorrelation and partial autocorrelation of longitude. (a) Longitude autocorrelation function, (b) longitude partial autocorrelation function.

Figure 2. Autocorrelation and partial autocorrelation of latitude. (a) Latitude autocorrelation function, (b) latitude partial autocorrelation function.

Figure 3. Google Earth motion track map longitude and latitude information. (a) Google Earth motion map. (b) Longitude and latitude coordinate information.

Figure 4. Comparison of measured data and fitted data of GPS. (a) The measured data and fitted data of longitude. (b) The measured data and fitted data of latitude.

Figure 5. The error sequence of GPS. (a) The error sequence of longitude. (b) The error sequence of latitude.

Figure 6. Autocorrelation function of residual sequence. (a) Autocorrelation function of longitude residual sequence. (b) Autocorrelation function of latitude residual sequence.

Figure 7. The error filtering results of longitude and latitude. (a) The error of longitude. (b) The error of latitude.

Figure 8. The longitude and latitude error curve of each model. (a) The error of longitude. (b) The error of latitude.

Table 1

The correlation function characteristics of the model.

Model	Autocorrelation Function	Partial Autocorrelation Function
MA	Truncation	Trailing
AR	Trailing	Truncation
ARMA	Trailing	Trailing

Table 2

Rectangular coordinates of a Gaussian projection plane.

Point Number	1	2	3	4	5	6
X-coordinates	322,186.148	322,185.5618	322,159.4055	322,159.4055	322,158.8192	322,183.8031
Y-coordinates	3,801,532.042	3,801,501.221	3,801,470.887	3,801,470.887	3,801,440.067	3,801,408.76
Heights	395	391	388	385	383	385
Point number	7	8	9	10	11	12
X-coordinates	322,183.2169	322,208.2011	322,233.1854	322,233.1854	322,233.1854	322,232.5994
Y-coordinates	3,801,377.94	3,801,346.633	3,801,315.327	3,801,315.327	3,801,315.327	3,801,284.506
Heights	386	384	382	384	384	385
Point number	13	14	15	16	17	18
X-coordinates	322,232.0134	322,231.4273	322,256.4121	322,255.8262	322,255.2402	322,255.2402
Y-coordinates	3,801,253.686	3,801,222.866	3,801,191.559	3,801,160.739	3,801,129.918	3,801,129.918
Heights	386	385	384	381	380	380
Point number	19	20	21	22	23	24
X-coordinates	322,254.6543	322,254.6543	322,228.4972	322,202.34	322,480.8022	322,480.8036
Y-coordinates	3,801,099.098	3,801,099.098	3,801,068.764	3,801,038.43	3,801,022.048	3,801,022.039
Heights	381	381	384	386	384	385

Table 3

Partial correlation function of longitude and latitude.

Order Number	Partial Correlation Function of Longitude	Partial Correlation Function of Latitude
1	−0.2150	−0.2805
2	−0.7232	−0.4044
3	−2.5410	0.6886
4	1.5189	−1.9369
5	0.5226	2.1126
6	0.6765	0.6765
7	0.4181	0.4181
8	0.2553	0.2553
9	0.1881	0.1745
10	0.1089	0.1090
11	0.1063	0.1012
12	−0.1034	0.0989
13	−0.0217	0.0976
14	0.0344	0.09675
15	0.10384	0.09243
16	−0.0531	0.0881

Table 4

The data of autoregressive coefficient and noise variance.

	Longitude		Latitude
Order Number	$\mathbf{Autoregressive}\mathbf{Coefficient}{\mathbf{\phi }}_{\mathbf{i}}$	$\mathbf{Noise}\mathbf{Variance}{\mathbf{\theta }}^{\mathbf{2}}$	$\mathbf{Autoregressive}\mathbf{Coefficient}{\mathbf{\phi }}_{\mathbf{i}}$	$\mathbf{Noise}\mathbf{Variance}{\mathbf{\theta }}^{\mathbf{2}}$
1	0.5750	6.07993 × 10−6	0.9146	8.0237 × 10−6
2	0.2136		0.0211
3	0.0846		0.0104
4	0.0722		0.0123
5	−0.0485		−0.0035
6	0.0074		−0.4411
7	0.0276		0.3617
8	0.0030		0.0561
9	0.0661		0.0032
10	0.0127		0.0064

Table 5

The original error sequence and filtered error sequence of longitude and latitude.

	Absolute Value of Maximum Error/m	Mean Square Error/m	Mean Absolute Error/m
Longitude original error sequence	1.8125	0.266574	0.179414
Longitude filtered error sequence	1.2597	0.242909	0.169197
Latitude original error sequence	1.594	0.206185	0.123304
Latitude filtered error sequence	1	0.185163	0.116725

Table 6

The longitude errors of different models.

	Absolute Value of Maximum Error/m	Mean Square Error/m	Mean Absolute Error/m
Original error	0.9474	0.209081	0.278591
AR (6)	0.5683	0.137044	0.177516
AR (7)	0.5633	0.137429	0.177733
AR (8)	0.5584	0.137194	0.177533
AR (9)	0.5614	0.137537	0.178879
AR (10)	0.5526	0.130864	0.176819
AR (11)	0.5605	0.137209	0.180477
AR (12)	0.5692	0.135520	0.179191
AR (13)	0.5528	0.135469	0.191706
AR (14)	0.5597	0.135521	0.177114
AR (15)	0.5628	0.137509	0.182409
AR (16)	0.5658	0.138854	0.186109

Table 7

The latitude errors of different models.

	Absolute Value of Maximum Error/m	Mean Square Error/m	Mean Absolute Error/m
Original error	0.9248	0.209406	0.270833
AR (6)	0.4301	0.127979	0.213698
AR (7)	0.4346	0.128573	0.212659
AR (8)	0.4379	0.129551	0.211016
AR (9)	0.4392	0.130159	0.211360
AR (10)	0.4258	0.126727	0.207835
AR (11)	0.4395	0.131173	0.209522
AR (12)	0.4322	0.128886	0.210099
AR (13)	0.4256	0.127212	0.209449
AR (14)	0.4416	0.127643	0.210880
AR (15)	0.4360	0.130439	0.212706
AR (16)	0.4546	0.133593	0.214826

References

1. Zhang, X.; Li, Z.; Wang, Z. Application of GPS positioning technology in different fields. Surv. Mapp. Inf. Eng.; 2001; 2, pp. 10-13.

2. Liu, X. Application of GPS survey technology in engineering survey. Sci. Technol. Innov. Appl.; 2020; 20, pp. 142-144.

3. Feng, Y.; Zhang, C. On the measurement principle and error analysis of GPS. J. Hebei Inst. Archit. Constr. Eng.; 2011; 29, pp. 94–96, 99.

4. Yang, B.; Yang, E. A Survey on Radio Frequency based Precise Localisation Technology for UAV in GPS-denied Environment. J. Intell. Robot. Syst.; 2021; 103, pp. 1-30. [DOI: https://dx.doi.org/10.1007/s10846-021-01500-4]

5. Ray, J.D.; Vijayan, M.S.M.; Kumar, A. Noise characteristics of GPS time series and their influence on velocity uncertainties. J. Earth Syst. Sci.; 2019; 128, pp. 1-14. [DOI: https://dx.doi.org/10.1007/s12040-019-1179-5]

6. Lei, K.; Ma, F.; Chen, B.; Luo, Y.; Cui, W.; Zhou, Y.; Liu, H.; Sha, T. Three-Dimensional Surface Deformation Characteristics Based on Time Series InSAR and GPS Technologies in Beijing, China. Remote Sens.; 2021; 13, 3964. [DOI: https://dx.doi.org/10.3390/rs13193964]

7. Wang, Y.; Zhao, J.; Fu, L.; Sun, X.; Chen, L. GPS positioning error analysis based on time series. Mod. Sci. Instrum.; 2012; 2, pp. 45-49.

8. Cardoso-Pereira, I.; Borges, J.B.; Barros, P.H.; Loureiro, A.F.; Rosso, O.A.; Ramos, H.S. Leveraging the self-transition probability of ordinal patterns transition network for transportation mode identification based on GPS data. Nonlinear Dyn.; 2022; 107, pp. 889-908. [DOI: https://dx.doi.org/10.1007/s11071-021-07059-x]

9. Dai, X.; Jiang, S.; Wu, C.; Jia, H. GPS time series analysis of Chengdu station based on the Lomb-Scargle algorithm. Arab. J. Geosci.; 2021; 14, pp. 1-8. [DOI: https://dx.doi.org/10.1007/s12517-021-08078-z]

10. Liu, Z.; Ma, X.; Giannakis, G.B. Space-time coding and Kalman filtering for time-selective fading channels. IEEE Trans. Commun.; 2002; 50, pp. 183-186.

11. Ho, T.J.; Chen, B.S. Tracking of dispersive DS-CDMA channels: An AR-embedded modified interacting multiple-model approach. IEEE Trans. Wirel. Commun.; 2007; 6, pp. 166-174. [DOI: https://dx.doi.org/10.1109/TWC.2007.05001]

12. Komninakis, C.; Fragouli, C.; Sayed, A.H.; Wesel, R.D. Multi-input multi-output fading channel tracking and equalization using Kalman estimation. IEEE Trans. Signal Process.; 2002; 50, pp. 1065-1076. [DOI: https://dx.doi.org/10.1109/78.995063]

13. Rong, W.; Si, C.; Weifeng, H.; Chun, D. Simulation analysis of satellite navigation and positioning error. Comput. Meas. Control; 2018; 26, pp. 264-268.

14. He, Y. Time series analysis modeling example. China New Commun.; 2020; 22, pp. 138-140.

15. Su, W.; Ran, S.; Liu, F.; Zhao, J. Multidimensional correlation time series modeling method based on correlation variable screening partial least squares regression. Inf. Control; 2021; 50, pp. 395-402.

16. Guo, Y. Accuracy analysis of Beidou GPS Integrated Navigation. Electron. Technol.; 2020; 49, pp. 58-59.

17. Feng, S.; Liu, X.; Wang, Y. Analysis of GNSS location time series by least square fitting. Surv. Mapp. Sci.; 2015; 40, pp. 157-160.

18. Specht, M. Statistical Distribution Analysis of Navigation Positioning System Errors-Issue of the Empirical Sample Size. Sensors; 2020; 20, 7144. [DOI: https://dx.doi.org/10.3390/s20247144]

19. Xu, W.; Si, Y. Analysis of noise characteristics of GPS time series—Taking Nepal and southern Tibet as an example. Surv. Mapp. Eng.; 2020; 29, pp. 23-30.

20. Zhang, Y.; Jiang, X.; Yang, W.; Liu, J.; Wei, Y.; Zhou, X. Research on joint positioning technology of GNSS system based on Kalman filter. Electron. Meas. Technol.; 2021; 44, pp. 109-113.

21. Peng, Y.; Chen, Y. Beidou GPS joint positioning algorithm based on Improved Kalman filter. Comput. Technol. Autom.; 2018; 37, pp. 59-62.

22. Zhang, G.; Han, S.; Ye, J.; Hao, R.; Zhang, J.; Li, X.; Jia, K. A method for precisely predicting satellite clock bias based on robust fitting of ARMA models. GPS Solut.; 2021; 26, 3. [DOI: https://dx.doi.org/10.1007/s10291-021-01182-3]

23. Meng, X.; Wang, Y. An improved adaptive extended Kalman filter algorithm for SINS/GPS compact integrated navigation system. J. Beijing Univ. Technol.; 2018; 38, pp. 625–630, 636.

24. Tian, J.; Liu, Y.; Wu, Y.; Lu, Y.; He, Z. GPS/INS Fusion Algorithm Based on Variational Bayesian Adaptive Kalman Filter. Int. Core J. Eng.; 2021; 7, pp. 26-32.

25. Bhattacharyya, S.; Mute, D. Kalman Filter-Based RAIM for Reliable Aircraft Positioning with GPS and NavICConstellations. Sensors; 2020; 20, 6606. [DOI: https://dx.doi.org/10.3390/s20226606]

26. Wang, Y.; Wang, H. Simulation Research on Kalman filtering of SINS/GPS integrated navigation system. J. Shenyang Univ. Technol.; 2019; 38, pp. 58-63.

27. Cao, S. Research on Integrated Navigation System and Filtering Algorithm Based on SINS/GPS; South China University of Technology: Guangzhou, China, 2016.

28. Specht, M. Consistency analysis of global positioning system position errors with typical statistical distributions. J. Navig.; 2021; 74, pp. 1201-1218. [DOI: https://dx.doi.org/10.1017/S0373463321000485]

29. Swaroop Kumar, N.; Ramesh, K.S. Analysis of carrier phase ambiguities to derive GPSTEC for precise positioning of GPS receivers. J. Crit. Rev.; 2019; 6, pp. 483-486.

30. Saracoglu, A.; Sanli, D.U. Accuracy of GPS positioning concerning Köppen-Geiger climate classification. Measurement; 2021; 181, 109629. [DOI: https://dx.doi.org/10.1016/j.measurement.2021.109629]