An Improved Divergence-Free Hatch Filter Algorithm Toward Sub-Meter Train Positioning With GNSS Single- Frequency Observations Only

Train positioning is the core function in the application of Global Navigation Satellite Systems (GNSS) in Railway Transportation. However, the use of the differential GPS (DGPS) along the Qinghai-Tibet Railway is expensive and difficult to maintain. Thus, a novel single-frequency algorithm based on the divergence-free Hatch filter is proposed, and no real-time augmentation correction input is required. The classical Hatch filter is severely affected by the divergence problem due to the ionospheric variation. In our algorithm, a novel decomposition-ensemble model is proposed for denoising and modeling the ionospheric variation, where the Variational Mode Decomposition (VMD) method is applied. With the aid of a sliding ionospheric variation fitting window, the divergence-free Hatch filter is constructed. The entire method is a so-called self-modeling method, but more efficient than recent studies. Besides, the Kalman filter is used for keeping continuous positioning accuracy. Finally, a static experiment in Tibet and a kinematic field test on the Qinghai-Tibet Railway is performed. In the ionospheric variation calculation-experiment, the experimental results show that the sliding window of our method can be shortened to 5 minutes with the data of 1s sampling rate, which basically meets the requirements of train positioning. In terms of train positioning accuracy, only the horizontal accuracy is concerned. In the static experiment, our method satisfies the accuracy requirements of the sub-meter level with a Root Mean Square Error (RMSE) value of better than 0.5m. In the kinematic test, the accuracy of our method is basically at the sub-meter level, with an RMSE value of approximately 0.6m.


I. INTRODUCTION
Global Navigation Satellite Systems (GNSS) have been diffusely applied in the domain of railway transportation, such as train control systems, railway fleet management [1], [2]. Train positioning, as precisely as possible, is the core requirement for the above GNSS applications. It should be noted that in train positioning, horizontal accuracy is mainly concerned [3]. At present, the differential GPS (DGPS) technology has been widely used in train positioning [4]- [7]. However, train positioning with the DGPS requires the The associate editor coordinating the review of this manuscript and approving it for publication was Zheng H. Zhu . support of a customized data network and a large number of base stations along the track. Therefore, the DGPS application leads to high construction costs and requires a lot of maintenance work, which is challenging to implement in harsh environments, such as the area along the Qinghai-Tibet Railway. All in all, a low-cost non-differential GNSS positioning algorithm is urgently required for train positioning. Appreciably, it is of considerable significance to attain sub-meter positioning through an inexpensive GNSS single-frequency receiver and a relatively simple algorithm.
In general, GNSS provides two kinds of fundamental observations, called pseudo-range observations and VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ carrier-phase observations. The pseudo-range observation is the pseudo-distance between a satellite and GNSS receiver, which is biased and coarse-ranged with a relatively large noise. Thus, its accuracy is reduced. However, the accuracy of the carrier-phase observation can attain the millimeter level. But it needs to fix ambiguities [8], [9]. Previous studies have demonstrated that by smoothing the pseudo-range observations with carrier-phase observations, it can lower the pseudo-range noise and enhance positioning accuracy [10]. Indeed, the phase-smoothed pseudo-range algorithm combines the advantage of these two observations, whereas significantly improving the accuracy of pseudo-range observations without introducing ambiguity. Among the existing smoothing algorithms, the Hatch filter algorithm is a classical method [11]. Yet the accuracy of the Hatch filter is seriously affected by the divergence in smoothed pseudo-range, which is caused by ionospheric variation. In view of the shortcomings of the Hatch filter, scholars have proposed several improvement strategies. An important strategy is to adjust the width of the smoothing window adaptively. For example, Park and Kee [12] introduced the Klobuchar ionospheric model to the calculation of the smoothing window. Zhang and Huang [13] applied external differential information from the satellite-based augmentation system (SBAS) rather than the Klobuchar model. Zhenggang et al. [14] comprehensively used satellite elevation angle, external augmentation information, and the distance change from the user to the base station to calculate the width of the smoothing window for a ground-based augmentation system (GBAS). Geng et al. [15] proposed the Three-Thresholds and Single-Difference (TT-SD) Hatch filter to achieve sub-meter positioning with Android devices. Zhou and Li [16] suggested using the Doppler phase difference observations rather than the carrier-phase observations to estimate the window width in kinematic positioning. Kim et al. [17] calculated the optimal window based on the ionospheric gradient estimation. However, the above approaches not only require frequent adjustment of the smoothing window, but the ionospheric bias is still included in the Hatch filter inevitably. Another strategy is to use external differential information for ionospheric correction to avoid divergence problem [18]- [22]. However, this strategy is challenging to implement under harsh conditions due to the reasons mentioned at the beginning of this paper, and the related correction algorithms are cumbersome. The authors are also inclined to solve this problem without any real-time augmentation correction input to reduce the cost and difficulty of implementation. Thus, this strategy is no longer discussed in detail. In fact, this problem can be solved easily with the GNSS dual-frequency receiver. Because the ionospheric variation can be easily calculated by dualfrequency observations, and the divergence-free Hatch filter model can be easily constructed with the calculated rate of ionospheric variation [23]- [25]. Of course, the above methods applicable to dual-frequency users don't apply to singlefrequency users.
Thus, Zhang et al. [26], Zhengsheng et al. [27] recently proposed an improved Hatch filter suitable for single-frequency users and didn't need any real-time augmentation correction input. In their studies, the concept of self-modeling was first proposed, that is, without any external information, only the fundamental single-frequency observations were used to model the ionospheric variation rate. Then a divergence-free Hatch filter model is built on this basis. However, there are still several shortcomings in this method. In particular, only a rough polynomial fitting method was used in the process of modeling and denoising the ionospheric variation rate. Thus, the width of the sliding ionospheric variation fitting window is long, and a 20-minute initialization time is required in real-time positioning to achieve the desired accuracy, which is unbearable for train positioning. Besides, this study used the ordinary least square method for parameter estimation, which makes it challenging to keep continuous positioning at the desired accuracy when the data quality is poor [15].
In order to improve the above method, this paper proposed an improved algorithm with a novel model to denoise and model the ionospheric variation, which is a type of decomposition-ensemble model [28]. Besides, the Variational Mode Decomposition (VMD) method [29] was introduced to construct this model, which has many advantages over empirical mode decomposition (EMD) or wavelets. The entire method is still a self-modeling method but more efficient. After the above improvements, the sliding ionospheric variation fitting window width can be shortened to 5 minutes, which basically meets the requirements of train positioning. Meanwhile, the Kalman filter model was introduced to maintain the continuity of positioning accuracy. Finally, a static experiment in Tibet and a kinematic field test on the Qinghai-Tibet Railway will be performed to verify the accuracy of our improved method and the effectiveness of the novel model to denoise and model the ionospheric variation.
The remainder of this paper is constituted as follows: Firstly, the classical Hatch filter with concise error analysis is reviewed in brief. The divergence-free Hatch filter model is discussed as well. Secondly, an improved algorithm is proposed. In this chapter, some related concepts and knowledge, i.e., the VMD method, the sliding window method, the self-modeling method are reviewed at first. Then, the improved model for denoising and modeling the ionospheric variation is discussed in detail. Next, the Kalman filter model for this algorithm is introduced. At the end of the chapter, we conclude with a summary of the improved algorithm. Thirdly, experimental results with real data are illustrated. Lastly, concluding remarks are presented at the end.

A. THE CLASSICAL HATCH FILTER
The Hatch filter is a simple recursive filter using the current observation and the former estimate without any extra sensors or dynamic models [30]. Thus, it can be applied to the realtime positioning with an inexpensive GNSS single-frequency receiver [22]. In the following, we will give a brief derivation of its mathematical model.
Ignoring multipath, tropospheric delay, and clock bias of the satellite and users' receiver, the GNSS observation equation at the kth epoch is expressed as: In the above equation, L k expresses the carrier-phase observation (denoted by distance). P k expresses the pseudo-distance observation. ρ k expresses the geometric distance between the satellite and the receiver. I k expresses the ionospheric delay (denoted by distance). N expresses the ambiguity (denoted by distance), which includes the hardware delay. ε L and ε ρ express the noise of observations. First of all, it should be noted that ambiguity and the ionospheric delay can be regarded as a constant within a period (k epochs) in the classical Hatch filter. This is because the ionospheric term varies very slowly with time [31].
The difference between the two observation formulas in equation (1) can be obtained as: Then a new definition is proposed as: In the above equation, A expresses the sum of ambiguity and twice the ionospheric delay in the starting epoch. According to the previous assumptions, ambiguity and the ionospheric delay can be regarded as unchanged. Thus, A can be estimated by simple averaging. Define A k as the average of A in k epochs, which can be expressed as: On the basis of equation (2) and (4), the smoothed pseudorange on epoch k can be obtained as: Thus, a recursive formula can be obtained simply by subtracting the adjacent two epoch elements, which can be expressed as: Substituting equation (4) into equation (6), the mathematical formula of the classical Hatch filter can be acquired: It should be noted that when cycle slips occur, the previous basic assumption is no longer valid and the Hatch filter must be reset.

B. HATCH FILTER ERROR CAUSED BY IONOSPHERIC VARIATION
According to the analysis in the previous section, the Hatch filter is a very concise method that can be easily implemented under general conditions. However, due to the influence of the ionospheric propagation delay, there will be a divergence problem when the ionosphere changes violently or after a long period of smoothing [22]. Thus, the influence of the ionosphere will be discussed in detail in this section.
In fact, unlike the previous assumption, the ionospheric delay isn't constant. Define I i,1 as the ionospheric delay difference between epoch i and epoch 1. The equation (4) can be converted into the following equation: According to the equation (8), There is an inherent bias between the estimated value A k and the actual value A. The ionospheric variation between epochs which take part in the smoothing leads to this bias. Thus, a divergence-free Hatch filter can be easily established on the basis of the correction of ionospheric variation, which can be expressed as: Now, the crux of the above problem is the calculation of the ionospheric variation rate.

C. CALCULATION OF THE RATE OF IONOSPHERIC VARIATION WITH DUAL-FREQUENCY OBSERVATIONS
As mentioned earlier, the ionospheric variation rate can be estimated with high accuracy by using dual-frequency carrier-phase observations.
Define L 1 and L 2 as the carrier-phase observations on two carrier frequencies (represented by distance). According to the equation (1), the difference between L 1 and L 2 can be obtained as: In the above equation, L k expresses the carrier-phase observation (represented by distance). b 1r , b 2r , b 1s , b 2s express the hardware delay on different frequencies of receiver/satellite, respectively. c TEC expresses the total electron contents (TEC) on the signal propagation route. f 1 and f 2 express the frequencies of the two carries. N 1 and N 2 express the ambiguity without hardware delay (represented by distance).
Then several new definitions are given as follows: In the above equation, I L 1 , I L 2 represent the ionospheric delay on different carriers. Then the equation (10) can be simplified as: (12) can be simplified as: In addition, the inter-frequency carrier-phase difference can be expressed as: B can be considered as a constant within a period (i.e., within the visible range of the satellite) [26], [27]. Then, both sides of the equation (14) are multiplied by the same coefficient. A new equation can be obtained as: For the GPS L1 signal: In the above equation, I L 1 ,k expresses the difference of ionospheric delay between the epochs. The difference between the carrier-phase observations on the two frequencies can eliminate the majority of the errors. Thus, the rate of ionospheric variation can be accurately determined through equation (16) with high sampling rate data, which can be used as the reference values for subsequent experiments [22], [26].

III. THE IMPROVED ALGORITHM
In order to help readers to understand the improved algorithm in this paper, this section will first briefly review some related concepts and knowledge, i.e., the VMD method, the sliding window method, and the self-modeling. Next, we will introduce the improvement points of the proposed algorithm in detail and, finally, conclude with the summary of the algorithm. The relevant experimental results and analysis will be given in Section IV.

A. VARIATIONAL MODE DECOMPOSITION
The VMD method is a novel signal decomposition method that decomposes a non-linear and non-stationary signal into discrete numbers of modes [32]. In essence, this method is substantially related to the generalized Wiener filter [29].
Compared with the widely used EMD method and wavelets, it has the following advantages: 1) Its theoretical basis is complete and has been strictly mathematically deduced. 2) It overcomes the mode mixing problem of the EMD method.
3) It overcomes the drawbacks of the end effect in the EMD method. 4) It has better noise and sampling robustness. 5) It overcomes the drawbacks of hard band-limits in wavelet methods. The procedure of using this method are mainly introduced below, and its detailed theory can be found in reference [29].
In the VMD method, the Intrinsic Mode Functions (IMF) are redefined as the amplitude-modulated-frequencymodulated (AM-FM) signals, which can be expressed as: In the above equation, the phase k is a monotone increasing function, and the envelope A k is non-negative. Both the envelope A k and the instantaneous frequency k change more gently than the phase k . The VMD method is aimed to decompose a complex signal into an ensemble of modes (IMFs) which are mostly compact around the corresponding center pulsations and to make the sum of the estimated bandwidth minimized. We can convert this problem into a constrained variational problem as follows: In the above equation, the u k (k = 1, 2, . . . , K ) and ω k (k = 1, 2, . . . , K ) are the set of modes and their center pulsations. The f expresses the original signal. In order to resolve the problem shown in equation (18), a quadratic penalty term α and Lagrange multiplier λ are applied to reconstruct the problem. In essence, this is a generalized Lagrange method, which can be represented as: In order to resolve the above non-constrained variational problem, the alternate direction method of multipliers (ADMM) is applied to the VMD method. The steps of the ADMM algorithm for VMD are as follows: initialize Dual ascent: ≤ ε output K modes and their center frequencies However, there are still two sub-problems in the above which need to be explained in detail.
First, the update process of the modes u k (k = 1, 2, . . . , K ) in the above can be rewritten as: Equation (20) describes a type of quadratic optimization problem. It is solved by the Parseval/Plancherel Fourier isometry under the L2 norm as follows: Second, the update process of the center frequencies ω k (k = 1, 2, . . . , K ) in the above steps can be rewritten as the following relevant problem: The optimization of this problem will be performed in the Fourier domain by the following equations (23) and (24): It should be noted that the estimated ω k is located at the power spectrum of the relevant mode, right in the center of gravity.

B. SELF-MODELING OF IONOSPHERIC VARIATION WITH SINGLE-FREQUENCY OBSERVATIONS
As discussed in Section II, the divergence problem of the classical Hatch filter is mainly caused by the ionospheric variation. The key to the divergence-free algorithm is modeling and denoising the ionospheric variation rate. For single-frequency users, this is challenging without any external information. As mentioned in the Introduction, the concept of self-modeling has recently been proposed to model and denoise ionospheric variation using only single-frequency observations. It was used to build an improved divergence-free Hatch filter suitable for single-frequency users without any real-time augmentation correction input [26], [27]. This significant concept will be reviewed in detail below.
According to equation (2), a new definition is proposed as: In the above equation, ε ρ k expresses the noise in pseudo-range observations. N expresses the ambiguity, which includes hardware delay (represented by distance). As discussed in Section II, the hardware delay can be considered as a constant. Therefore, when no cycle slips occur, the ionospheric variation is straightly synchronous with the variation of the I N . Consequently, we can calculate the ionospheric variation roughly by equation (25).
Obviously, under the influence of the pseudo-range noise, the ionospheric variation calculated by this method is not accurate. To solve this problem, Zhang et al. used a secondorder least-square polynomial fit method to estimate it with the observations of multiple epochs. Meanwhile, a sliding window was used to select the observations for estimation, that is, the sliding ionospheric variation fitting window. The entire method is named as the sliding window fitting method [26], [27]. It should be noted that the extrapolated value of the fitted polynomial is used when a cycle slip occurs.
The schematic diagram for the sliding window is illustrated in Figure 1.
As shown in Fig. 1, each small cell denotes the data of an epoch. Meanwhile, the cell filled with green represents the data of the current epoch. The sliding window is a data block for fitting, which consists of a fixed number of small cells. In data processing, the window keeps moving forward to process the entire data. For post-processing positioning, the current epoch is at the center of the window, and there is overlap between windows, which is called the segmented sliding window with overlap. For real-time positioning, the current epoch is at the end of the window, called the epoch sliding window [26], [27].
According to the studies of Zhang et al., on the basis of the above method, the size of the sliding window should be at least 20 minutes to achieve positioning towards sub-meter accuracy at the sampling rate of 1s [26], [27]. As mentioned in the Introduction, it is unacceptable in real-time train positioning. Thus, in the next subsection, we will propose an improved model to denoise and model the ionospheric variation.

C. IMPROVED MODEL TO DENOISE AND MODEL THE IONOSPHERIC VARIATION
As mentioned above, an improved model is proposed in this subsection, which is a typical decomposition-ensemble model. This concept has recently been used to explain why multi-scale complexity decomposition can improve fitting and predictions [28]. The framework of the decompositionensemble model for fitting is shown in Figure 2.
As mentioned in the former subsection, a second-order least-square polynomial fit method was used to fit the ionospheric variation. Based on it, the VMD method whose superiority has been fully confirmed in existing studies is introduced to construct the improved model, which is a typical decomposition-ensemble model.
As shown in Fig.2, the improved model in this study is composed of three main components, i.e., multi-scale decomposition, individual fitting and ensemble fitting, which are described in detail as follows: 1) Multi-scale decomposition: In this study, the VMD method is employed to decompose the complex original data into several relatively simple modes, which is a multi-scale analysis technique. 2) Individual fitting: all modes obtained can be fitted by a specific fitting algorithm, which is the second-order least-square polynomial fit algorithm in this study. Several individual results are obtained. 3) Ensemble fitting: combine the individual results and reconstruct them into the ensemble result. The reason why the decomposition-ensemble model introduced above can improve the fitting accuracy is that it can mitigate modeling difficulty [33]. More specifically, the complex original data are converted into relatively simple modes through this model. Thus, the modeling complexity of original data can be reduced. According to the complex system theory, a system that has a lower complexity is more inclined to follow a deterministic process. Thus, it can be captured and modeled more easily [34]. As shown in Fig.2, it converts a difficult fitting task for original data with high complexity into several manageable fitting subtasks with low complexity. Thereby, the ensemble fitting result with higher accuracy can be obtained. On the other hand, the inherent driving factors in complex non-linear systems can be captured using the VMD method (in terms of IMF), which can make the obtained modes simpler and more meaningful (in terms of low complexity) [29], [35]. Besides, the VMD method is greatly related to Wiener filter denoising. It can perform a few denoising effects by setting a sizeable quadratic penalty term and shutting off the Lagrange multiplier.
All in all, this improved model is one of the significant improvement points compared to previous researches. It is used to model and denoise the ionospheric variation rate with the sliding window method introduced in the former subsection. The entire method is still a self-modeling method. Thus, the divergence-free Hatch filter suitable for single-frequency users can be constructed based on the ionospheric variation rate calculated above. After the above improvements, the width of the sliding ionospheric variation fitting window can be shortened to 5 minutes. The relevant experimental results and analysis are given in Section IV.

D. KALMAN FILTER
As mentioned above, in the process of smoothing pseudo-range with the proposed improved Hatch filter, when cycle slips occur, the Hatch filter must be reset. Thus, the continuity between epochs will be destroyed in terms of the position parameters. However, in the studies of Zhou et al., only the ordinary least-square method was used for parameter estimation. Therefore, it is challenging to guarantee continuous positioning accuracy in harsh environments or with imperfect data.
Thus, the Kalman filter is introduced for our algorithm. It is another improvement point compared to the studies of Zhou et al., and a cycle slip detection should also be considered in our algorithm. Refer to the reference [22], the timedifference of the code-minus-carrier can be used to detect the cycle slip. In particular, since the noise of pseudo-range has been smoothed, the combination of the filtered pseudo-range and carrier-phase is an excellent cycle slip indicator [22]. It should be noted that the filter must be reset once the cycle slip occurs.
Next, we will introduce the Kalman filter model used in this paper. Since this method is commonly used and classic, in order to make this paper concise, we only briefly introduce the key part.
Considering the safety and comfort of the train, the speed of the train changes relatively gently with a relatively small acceleration [36]. Therefore, a relatively concise and easy-to-implement eight-parameter model [37], [38] will be applied to our algorithm. The recursive state equation is: In the above equation, k expresses the kth observations. X expresses the state vector of the user. φ expresses the state transition matrix between epochs. w expresses the processing noise vector. Then, define T as the sampling interval and Q as the covariance matrix of w. The vector X is described as follows: In the above equation, the vector X consists of position components, velocities, clock bias δt and clock frequency drift δf . The transition matrix φ is described as follows: The processing noise matrix Q can be described as follows (29), as shown at the bottom of the next page, In the above equation, the variables Sv, S t , S f represent the spectral density of velocities, clock bias, and clock frequency drift.

E. SUMMARY OF THE IMPROVED ALGORITHM
We will summarize our improved algorithm in this subsection. The flowchart of it is shown in Figure 3. In detail, the processing procedures of the improved algorithm includes: 1) Input the single-frequency original observations data. 2) As mentioned in the C part of Section III, the sliding ionospheric variation fitting window is established for the self-modeling of ionospheric variation. In order to accomplish it, an improved model using the VMD method and the second-order least-square polynomial method is constructed, which is a typical decomposition-ensemble model for fitting. The entire method is still a self-modeling method. Thus, the ionospheric variation rate is obtained. 3) As mentioned in the B part of Section II and the D part of Section II, the divergence-free Hatch filter based on the above results is constructed to smooth the pseudorange. And the filter must be reset when a cycle slip occurs. 4) As mentioned in the D part of Section III, the Kalman filter method is introduced for parameter estimation. 5) The above steps can be cycled for calculation.
According to the reference [26], [27], several pre-stored precision correction products for other errors are used for positioning solutions, such as the prediction products for orbits, clock bias, and ionospheric delay, etc. It should be noted that the above precision products are stored in advance rather than real-time input precision correction information. There is no need to set up additional channels to receive them. Therefore, the application of these precision products is in line with our requirements. Thus, the Ultra-Rapid orbits and clock bias products, the VMF1 forecast products [39], and the global ionosphere model prediction products (e.g.C2PG products) are used for error corrections in this paper.

IV. EXPERIMENT AND ANALYSIS
In this section, the data of a static experiment in Tibet and a kinematic field test on the Qinghai-Tibet Railway are processed and analyzed. Thus, the effectiveness of the improved model for denoising and modeling the ionospheric variation is verified. And the positioning accuracy of the proposed improved algorithm is demonstrated.
The static experiment is a technical verification experiment conducted in advance. Since this paper mainly focuses on the area along the Qinghai-Tibet Railway, the L1 band data of a 1s sampling rate from the LHAZ station in Tibet are used for the static experiment in this paper. We take a four-hour period of data from 12:00 to 16:00 UTC with a satellite cutoff altitude of 15 degrees as an example for analysis.
The kinematic experiment was carried out on January 2, 2019. A test train for data collection started from NaQu Railway Station and traveled northward along the Qinghai-Tibet Railway, shown in Figure 4. A one-hour period of data from 7:00 to 8:00 UTC with a 1s sampling rate was collected for processing and analysis. The obtained data are used for the simulated real-time train positioning experiments.  As shown in Fig.4, the receiver antenna was fixed on top of the test train. And the CHC N72 GNSS receiver used in our experiment was located in the compartment of the test train, which is shown in Figure 5. To obtain the precise track of the test train as the reference true value of the kinematic experiment, we got it based on the digital tracking map and information from several location sensors provided by the relevant departments.
The reference motion trajectory of the test train is shown in Figure 6.
As mentioned in the Introduction, in train positioning, horizontal positioning accuracy is mainly concerned. Thus, the train positioning deviation is given under a Gauß-Krüger coordinate, denoted as δ H [3]. A schematic diagram of this deviation is illustrated in Figure 7.
In the remainder of this section, we will introduce experimental results and relevant analysis in three parts, respectively.

A. EXPERIMENTS ON IONOSPHERIC VARIATION
In this subsection, the results of experiments on ionospheric variation are shown and analyzed. As mentioned in Section III, we proposed a novel model to denoise and model the ionospheric variation. The improvement effect of the proposed model is demonstrated through the following experiments.
We use the four-hour period of data of the LHAZ station to estimate ionospheric variation by different approaches mentioned above. According to the dual-frequency calculation method in the C part of Section II, the ionospheric variation (calculated by Eq. (16)) of different satellites during the  experimental time is illustrated in Figure 8. It shows that the ionospheric variation between epochs is small, which also confirms with the previous conclusion.
As mentioned in the C part of Section II, the ionospheric variation estimated by the dual-frequency method can be used as the reference values of our experiments [22], [26].
Subsequently, the ionospheric variation is calculated by three different approaches, i.e., calculating raw results roughly by Eq. (25), the original self-modeling method [26], [27], and the improved method in this paper. As an example, we will show the result of the satellite G01, which is visible during the whole experimental period. Figures 9 and Figure 10 show the estimated residuals of the above three methods under the 5-minute and 20-minute sliding ionospheric variation fitting window, respectively. We also calculated the RMSE of the above three methods under the sliding ionospheric variation fitting window of various widths, as shown in Figure 11. Meanwhile, the values of them are given in Table 1.
According to the above results, we can find that the estimated residuals of the improved method with a 5-minute sliding ionospheric variation fitting window and the original self-modeling method with a 20-minute window are basically equivalent. That is to say, compared with the original method, the improved method can shorten the width of the sliding ionospheric variation fitting window to 5 minutes, which basically meets the requirements of train positioning. This shows that our method is very effective and practical.

B. STATIC EXPERIMENTS
The configurations of the static experiment have been introduced in detail at the beginning of this section. As mentioned above, we will focus on the horizontal accuracy of positioning and the accuracy in the north and east directions. We will use three methods for comparative experiments, i.e., positioning with raw observations, the classical Hatch filter, and the improved method with a 5-minute window. In order to carry out comparative experiments, it should be noted that the above three methods use the same error correction configurations (mentioned in the E part of Section III), except for pseudo-range observations. In other words, the above three methods all use pre-stored forecasting precision products to correct the main error terms.
Firstly, compared with the positioning results of raw observations, the positioning errors of the improved method are given. Figures 12 and Figure 13 show the positioning errors of the two methods in the east and north directions, respectively. Figure 14 shows the Horizontal positioning errors (i.e., the deviation illustrated in Fig.7) of the two methods varying with time. Figure 15 shows the Horizontal positioning errors under a Gauß-Krüger coordinate. Figure 16 shows the positioning errors of the classical Hatch filter, which has an obvious divergence problem.    In addition, we calculated the RMSE values of the positioning errors with three methods, which can be used to evaluate positioning accuracy. They are illustrated in Table 2.
According to the above results, we can find that the accuracy has been improved by our method compared to the  positioning with raw observations. The horizontal accuracy is better than 0.5m, which has significantly met the requirements of sub-meter positioning. And it shows that a 5-minute window is enough. Meanwhile, the control experiment shows significant noise effects. Of course, this is the conclusion obtained when the environment is relatively ideal, and the data quality is relatively good. After all, it is just a technically validating experiment. In addition, we can also find that the accuracy of the classical Hatch filter is significantly influenced by the divergence problem, which is caused by ionospheric variation (discussed in Section III).

C. KINEMATIC EXPERIMENTS
The configurations of the kinematic experiment have been introduced in detail at the beginning of this section. This is    Figure 18 show the positioning errors of the two methods in the east and north directions, respectively. Figure 19 shows the Horizontal positioning errors of the two methods varying with time. Figure 20 shows the positioning errors of the classical Hatch filter, which also has an obvious divergence problem.
Similarly, we calculated the RMSE values of the positioning errors with three methods, which can be used for evaluating positioning accuracy. They are illustrated in Table 3.
According to the above results, we can find that the accuracy has been improved by the improved method compared to the results with raw observations in train positioning. As shown in Fig.19, the horizontal accuracy can converge to the sub-meter level in a very short time, and the overall errors are less than 1m. By contrast, a lot of horizontal errors of raw observations exceed 1m due to the noise, and the control experiment shows significant noise effects. In addition,    we can also find that the accuracy of the classical Hatch filter is significantly influenced by the divergence problem. All in all, our algorithm basically meets the requirements of sub-meter positioning. However, when the environment is harsh, or the data quality is poor, there are still a small number of outliers inevitably. Of course, this can be easily solved by digital map matching, inertial navigation assistance, and other mature technologies, which is not the focus of this paper. Generally speaking, this method basically meets the requirements for GNSS in train positioning.

V. CONCLUSION
In the domain of train positioning, a low-cost non-differential GNSS positioning algorithm is urgently required. The improved Hatch filter algorithm using an inexpensive singlefrequency receiver is a meaningful idea. The Hatch filter has been one of the most frequently-used filters due to the simplicity of construction and its effectiveness in reducing noise. It combines the advantage of the pseudo-range and carrier-phase observations, whereas significantly improving the accuracy of pseudo-range observations without introducing ambiguity. However, the accuracy of the classical Hatch filter is severely affected by the divergence problem due to the ionospheric variation. In order to construct the divergence-free Hatch filter, the key is to calculate ionospheric variation, which is challenging for single-frequency users without any real-time augmentation correction input.
In this paper, we proposed a novel algorithm based on the divergence-free Hatch filter using single-frequency observations only, which does not need any real-time augmentation correction input. To be exact, this method was produced by further improvement on the basis of the original self-modeling method proposed in recent studies. Compared to the ordinary method, three innovations were included in this method: 1) A novel model with the sliding ionospheric variation fitting window for denoising and modeling the ionospheric variation, which is a type of decomposition-ensemble model for fitting. The entire method is still a self-modeling method, but more efficient than the ordinary method in recent studies. 2) The emerging VMD method for decomposition, which is significantly related to the generalized Wiener filter. It is also helpful for denoising. 3) The Kalman filter used for parameter estimation to keep continuous positioning at the desired accuracy.
In this paper, a static experiment in Tibet and a kinematic field test on the Qinghai-Tibet Railway was performed to verify the accuracy and the effectiveness to denoise and model the ionospheric variation of our algorithm. The experimental results show that the estimated residuals of our improved method with a 5-minute sliding ionospheric variation fitting window and the original self-modeling method proposed in recent studies with a 20-minute window are basically equivalent. It illustrated that the sliding window width could be shortened to 5 minutes with the data of 1s sampling rate, which basically meets the requirements of train positioning.
In the verification experiments of positioning accuracy, three methods were conducted for comparative experiments, i.e., positioning with raw observations, the classical Hatch filter, and our method with a 5-minute window. It should be noted that the above three methods use the same error correction configurations (mentioned in the E part of Section III), except for pseudo-range observations. In the static experiment, our method satisfies the accuracy requirements of the sub-meter level with an RMSE value of less than 0.5m. Meanwhile, the control experiment showed significant noise effects. However, this conclusion was obtained when the environment is relatively ideal, and the data quality is relatively good. After all, it was only a technically validating experiment. In addition, we found that the accuracy of the classical Hatch filter is greatly influenced by the divergence problem. In the kinematic experiment, the procedure of the comparative experiments was similar to the static experiment. And the accuracy of our method was basically at the sub-meter level with an RMSE value of approximately 0.6m. Similarly, the control experiment showed significant noise effects, and the classical Hatch filter showed a serious divergence problem. All in all, our method basically met the requirements for the application of GNSS in train positioning.