A Novel Dynamical Filter Based on Multi-Epochs Least-Squares to Integrate the Carrier Phase and Pseudorange Observation for GNSS Measurement

: The high noise of pseudorange and the ambiguity of carrier phase observation restrain the GNSS (Global Navigation Satellite System) application in military, industrial, and agricultural, to name a few. Thus, it is crucial for GNSS technology to integrate the pseudorange and carrier phase observations. However, the traditional method proposed by Hatch has obtained only a low convergence speed and precision. For higher convergence speed and precision of the smoothed pseudorange, aiming to improve positioning accuracy and expand the application of GNSS, we introduced a new method named MELS (Multi-Epochs Least-Squares) that considered the cross-correlation of the estimating parameters inspired by DELS (Double-Epochs Least-Square). In this study, the ionospheric delay was compensated, and so its impact was limited to the performance of the ﬁlters, and then exploited the various ﬁlters to integrate carrier phase observation and pseudorange. We compared the various types of Hatch’s ﬁlter and LS (Least-Square) methods using simulation datasets, which conﬁrmed that the types of LS method provided a smaller residual error and a faster convergence speed than Hatch’s method under various precisions of raw pseudorange. The experimental results from the measured GNSS data showed that LS methods provided better performance than Hatch’s methods at E and U directions and a lower accuracy at N direction. Nevertheless, the types of LS method and Hatch’s methods improved about 12% and 9–10% at the 3D direction, respectively, which illustrated the accumulating improvement at the enhanced directions was more than the decreased direction, proving that the types of LS method resulted to better performance than the Hatch’s ﬁlters. Additionally, the curve of residual and precision based on various LS methods illustrated that the MELS only provided a millimeter accuracy di ﬀ erence compared with DELS, which was proved by the simulated and measured GNSS datasets.


Introduction
The GNSS (Global Navigation Satellite System) technology has been widely used for precision agriculture, time synchronization and delivery, and deformation monitoring, to name a few, although the systems are still constructing and upgrading. As of February 2020, GPS (Global Positioning System) with 31 operational satellites are still undergoing modernization, which is maintained and GPS and did not provide references for alone-receiver. Most of the scholars have demonstrated their methodology with an indirect method in lights of positioning results because it is hard to provide true pseudorange, and thereby Qian et al. [24] simulated GNSS observations by adding diverse random error to the real satellite-to-receiver distance. They obtained a higher precision smoothed pseudorange than Hatch's smoother by employing the recursive LS, which is only used for an ionospheric removed GNSS. To counteract the influence of ionospheric delay and improve the precision of smoothed pseudorange, Chen [25] et al. exploited a similar solution to the recursive LS [24] to integrate the pseudorange and TDCP (Time Difference Carrier Phase), which divided GNSS measurements into dispersion and non-dispersion term. Then, they exposed the merit according to the residual error and positioning accuracy, which lacked analyzing convergence speed and precision, and application scope. The presented methodologies in the light of parameter estimation only provided a low filtering accuracy and convergence speed, and even some contribution lacked analyzing the two critical indicators in most cases.
To achieve higher integration accuracy and convergence speed, thereby improving positioning accuracy and expanding the scope of GNSS applications, this study introduced a recursive MELS (Multi-Epochs Least-Squares) data processing strategy to integrate the RP (Raw Pseudorange) and TDCP measurements. The developed algorithm has the advantage of requiring the same information compared to the existing integration methods, and it is easy to implement and graft to any projects with an enormous amount of measurements because of a small calculation. This study firstly discussed the details on the GNSS positioning principle and techniques to incorporate the observation, such as Hatch's filtering and the existing LS integration approach, and then presented an enhanced recursive MELS method for the fusing phase carrier and pseudorange observations. The effectiveness of the developed method was proved through the analysis of simulated dataset and field measurements and, finally, was discussed the merit of the presented method and its contribution to body knowledge of the carrier phase smoothing pseudorange technology.

Traditional Integration Method
This section introduces the fundamental knowledge to integrate pseudorange and phase carrier observation for multi-frequency GNSS, and then the traditional Hatch's filter in terms of its precision has been discussed.

Mathematical Models
For a ground-based alone-receiver r receiving signals from the observed satellite j, the observation of pseudorange P and carrier phase L in terms of distance for multi-frequencies GNSS can be expressed as, P j i,r = ρ j r + c · (dt j − dτ r ) + I j i,r + T j r + ε j p,i,r where the superscript j is the satellite index, while the subscript i and r refers to the index of signal frequency and receiver, respectively; the ρ represents geometric distance as a function of the receiver and satellite coordinates (m); the c is the speed of light in a vacuum (m/s); the dt and dτ are the clock offset of satellite and receiver, respectively (s); the I represents ionospheric delay based on the relationship with the frequency, where f is the frequency of the signal (Hz); the T indicates tropospheric delay (m); λ denotes the wavelength of the signal (m); N represents integer ambiguity (cycle); the ε is the residual error of carrier phase and pseudorange measurements (m). The instrument hardware delay is ignored in this study because of its small influence. Moreover, the PCO (Phase Center Offset) and PCV (Phase Center Variation) of receiver and satellites, phase wind-up correction, earth rotation correction, relativistic effect, tidal loading, and other modellable errors make an essential influence on the positioning performance of GNSS, all of which can be compensated by Remote Sens. 2020, 12,1762 4 of 20 the corresponding model for a more accurate result [26][27][28][29][30]. Besides, the precision products provided by IGS (International GNSS Service) organization can suppress satellite orbit error and clock offset, and the residual errors can be ignored in practical engineering applications.
When we applied the TDCP technology for integrating the carrier phase observation and pseudorange, some errors had a small influence on the filtering performance, such as the PCO and PCV of receiver and satellites, the phase wind-up correction, the earth rotation correction, the relativistic effect, tidal loadings, and other modellable errors. Ignoring the errors mentioned above because of their slowly changing speed, we illustrated a time difference form relying on Equations (1) and (2), and it is the following: where the t and t − 1 is the measuring time as the receiver is collecting GNSS signals; the prefix ∆ represents TD (Time Difference) symbol, with ∆P and ∆L (also called TDCP observation) serving the increment variation of the two typical measurements as time increases from epoch t − 1 to t, respectively. The TDCP measurements between two close epochs eliminate the integer ambiguities if there is no-existence or have repaired cycle slip [31], as well as suppressing satellite orbit error and clock offsets, and the modellable errors because of their slowly changing speed and small influence. On the other hand, the non-dispersive terms make the same influence on the two typical GNSS observations, thereby summarizing as part of the receiver-to-satellite distance. Thus, the errors make no influence on the final integration solution. However, the ionosphere delay is a significant error for all integrating methods due to its strong impact [32], which would decrease the precision of IP (Integrated Pseudorange) by accumulating the distance between carrier phase observation and pseudorange. The ionosphere variation in Equations (3) and (4) is an unignored value for integrating carrier phase and pseudorange, and we compensated for the error rely on the IF (Ionospheric-Free) combination, as well as its variation, for dual-frequency receivers. Then, Equations (5) and (6) illustrate simplified TD form of Equations (1) and (2) for an ionosphere-removed GNSS as following: where the subscript i f indicates ionosphere-free measurements, and the ∆ ρ denotes the satellite-to-receiver distance variation, considering the non-dispersion errors of tropospheric delay and the clock drift of receiver, which takes a value with ∆ ρ = ∆ρ − c · ∆dτ + ∆T.
Because we have compensated for the ionospheric delay error; thus, we no longer distinguished that the ionospheric delay was free/exist with removing the subscript, i f , for P and L in the following derivation.

The Smoother of Hatch
The Hatch's filter, as well as the expanded method, is the first type of phase carrier observation, smoothing pseudorange approach, which is well-known for its convenience and simplicity. Hatch's method relies on a quiet ionosphere with ∆I(t, t − 1) = 0, and the main idea is that the increment of pseudorange and carrier phase observation is equal to each other within a short period [7]. Hatch established an iterative form for the traditional method as follows, Remote Sens. 2020, 12, 1762 5 of 20 where theP is the smoothed pseudorange and k is the smoothing time constant, whose value is recommended to take a value with 20-100 for 1 s sampling interval. Because the convergence speed and accuracy are two main indexes to evaluate the performance of filters, the variance-covariance theorem provides the theoretical precision of the smoothed pseudorange of Hatch's smoother as follows, where the Q represents the variance of corresponding measurements decided by subscript. The initialization of the smoothed pseudorange is a critical problem to be processed for Hatch's smoother because the deviation of the initial value requires lots of epochs to approach the actual value. Besides, the smoother has suggested that it can obtain the initial smoothed pseudorange by the average of RP at a defined time t 0 , giving an equation as following, Next, the error propagation law provided the variance-covariance information for the initial value of the smoothed pseudorange according to Equation (10), taking a form as follows, Equations (7)-(11) is the iteration process of Hatch's smoother, which provides the value of IP and its variance information at any epoch for the subsequent accuracy analysis and comparison. Moreover, another urgent problem still to be processed is how to get an appropriate smoothing time constant because a small k is unable to shrink noise, and a significant k has accumulated much ionospheric variation error. Thus, the next section provides the theoretical accuracy of the Hatch method as changing the pseudorange accuracy, the time constant, and the filter time. To overcome the influence of the ionosphere delay, some scholars have applied the IF combination to the traditional method.

Least-Squares Fusion Algorithm
This section introduces the DELS (Double-Epochs LS) and analyzes its theoretical accuracy, and then the multi-epochs based on triple-epochs and quadra-epochs is also presented, which is defined as TELS (Triple-Epochs LS) and QELS (Quadra-Epochs LS), respectively.

Double-Epochs Least-Squares Fusion Algorithm
A new type of algorithm, named recursive LS parameter estimation methods, is also used for integrating the carrier phase observation and pseudorange, which has constructed a cost function to integrate typical measurements and requires prior accuracy information. The LS integrating method has designed the observation form for the raw measurements as Equation (12).
where the Y(t) expresses the vector-matrix of raw observations. Remote Sens. 2020, 12, 1762 6 of 20 Besides, according to the Equation (12), the Equation (13) has provided the variance-covariance matrix information of the virtual observation, Y(t), for the following theoretical precision derivation, where the Q Y(t) is the variance-covariance matrix of Y(t), relying on the prior information provided; theρ is the estimated pseudorange by the LS fusion algorithm; the t and t − 1 are the observation epoch of current and previous, respectively. Next, the LS parameter estimation method has constructed a virtual observation, relying on TDCP and RP observation by transformation matrices, T, and has deduced the corresponding variance-covariance matrix information according to the error propagation theorem, which is as follows: where the y(t) is the virtual observation and T is the transformation matrix. On combining the Equations (12)- (15), the Equation (16) has designed another form for the y(t), which is the core of the LS method for integrating carrier phase observations and pseudorange, which is expressed as follows, In the field of LS problem, y(t) delegates measurements; x(t) is the estimating parameter with a specific value ofρ(t), which is also called IP; H denotes the designed matrix with an invariant form, 1 1 T , at this question. Equations (12)- (16) are the fundamental equations of the LS integrating method, and then we can obtain the LS solution of the estimated value, based on the cost function, which is as follows, Equation (18) provides the theoretical variance-covariance of Equation (17) based on the error propagation theorem.
To achieve a rigorous parameter estimation method, the LS has considered the self-correlation among the raw observations and estimated parameters to construct the variance-covariance matrix for all of the epochs. According to the variance-covariance propagation theorem, Equation (19) gives the covariance information of the estimated value, x(t), and original measurement, Y(t), based on the relationship provided by Equation (17).
Remote Sens. 2020, 12, 1762 7 of 20 Equations (13), (18), and (19) clearly state that one of the matrix functions, F, can transform Q x(t−1) and Q x(t−1)Y(t−1) to Q Y(t) , which is the basis to deduce the theoretical convergence accuracy for LS integration method. This manuscript did not focus on the specific form of the transformation matrices or matrix functions and its equivalence and showed that any type of F is applicable as long as it possesses capability mentioned above.
The DELS method does not integrate the pseudorange and carrier phase observation at the first epoch. At the second epoch, the Q Y(t) is a diagonal matrix, and then the DELS filter can enter the loop by running Equations (12)- (20), which is the entire process of the DELS.
To obtain the theoretical convergence accuracy and make an in-depth analysis of the LS integration method, we assumed that it works as long as enough and becomes steady. It would establish the following relationship for the variance-covariance information of the current and previous epoch.
Based on the matrix functions provided by combining the Equations (12)- (22), the numerical accuracy convergence of the DELS method for integrating the carrier phase observation and pseudorange can be deduced. The GNSS observation is always limited, so we did not provide a specific form of convergence accuracy.

Multi-Epochs Least-Squares Fusion Algorithm
This section presents an enhanced recursive LS based on multi-epoch to integrate the carrier phase observation and pseudorange, where the TELS and QELS are presented as an example. The vectors mentioned in this section are all for multi-epoch, and we omitted the subscripts for the simplicity of formula expression.
Inspired by the original measure value expressed by Equation (12), the TELS method has enhanced the DELS by adding an epoch at one-time filtering for GNSS raw measurements, whose raw observation vector is expressed as Equation (23).
where the meaning of the symbol is the same as explained in Equation (12).
This manuscript enhanced the DELS by adding epoch at one-time filtering, which forms a new virtual observation as illustrated in Equation (25), whose variance-covariance matrix information has Remote Sens. 2020, 12, 1762 8 of 20 the same calculation method as that of Equation (15) but an unequal measure because the size and value of Equation (23) are both different from that of Equation (12).
Similar to Equation (16), Equation (26) has expressed the y(t), relying on the estimated values and observations of current and previous two epochs as reported by bringing Equations (23)-(25) together, which is the measurement equation of the recursive TELS integrating method and is as follows, In the process of solving the TELS problem, the meaning of the parameter is the same as in Equation (16) of the previous section, whereas the invariant design matrix provides a new value with Then, the LS method gives the numerical solution of x(t) as shown for Equation (17), which has the same form and different sizes with the TELS method. Besides, the error propagation law provides the variance-covariance information of the x(t) and Y(t), shown as Equations (18) and (19). The variance-covariance of the x(t) and Y(t) in the TELS method has the same form and different sizes with the DELS. From the Equations (23)- (26) and their related equations, the users can construct the stochastic model of Equation (23) for the next epoch, relying on Equations (18) and (19), if they have obtained the prior precision information of pseudorange and carrier phase observations. Similar to the DELS method, the summary of the above equations shows it can deduce the The TELS method does not integrate the pseudorange and carrier phase observation at the first epoch. At the second epoch, the Q Y(t) is a diagonal matrix, and the DELS is used to integrate the measurements. Then, the TELS filter can enter the loop by running Equations (23)- (27) and related equations, which is the entire process of the TELS.
To obtain the theoretical convergence accuracy and compare the performance of TELS with others, we established the following relationship for the variance-covariance information of the current and previous epochs, assuming that the filter works as long as enough and becomes steady.
Based on Equations (23)-(28), the users can obtain the theoretical numerical convergence accuracy of the TELS method for integrating the carrier phase observation and pseudorange.
The QELS has enhanced the DELS by adding two epochs at one-time filtering for GNSS raw measurements, whose raw observation vector is expressed in Equation (29) as follows, In this study, the specific result for QELS was not deduced, but this could be produced by deducing the numerical solution based on four or more epochs to estimate the IP, as well as its variance-covariance information. Furthermore, the comparison of the results between Equations (12)- (22) and Equations (23)- (28) proved that the DELS based on double-epochs was Remote Sens. 2020, 12, 1762 9 of 20 a particular case of the T/MELS, relying on triple-epochs, and the T/MELS was an expansion of the DELS.

The Comparison of Hatch's Filter and Least-Squares Methods
Because the existence of system errors causes challenges to obtaining the real pseudorange, which is essential information to evaluate the related techniques [24], a simulation dataset was used to compare the different methods for integrating the carrier phase observation and pseudorange. The emulator simulated different types of GNSS observations by adding various noise levels, which added 3 mm Gaussian white noise for carrier phase and series noise for RP to the satellite-to-receiver distance, relying on the final precise orbit of the satellite and the position of the receiver provided by IGS organization. The emulator outputted measurements at DOY (day of the year) 300 (27 October), 2018 for ALBH station, which was located in Canada. Based on the reference provided and satellite-to-receiver distance mentioned above, it took the E02 (the GALILEO navigation satellite with pseudo-random noise code of 2) as an example to test and verify the performance of Hatch's smoother, the DELS, and the M/TELS method. The subgraph A-F of Figure 1 shows the residual error of various methods as changing the precision of pseudorange accuracy, where the k represents the smoothing time constant of Hatch's smoother.
The residual error presented in subgraph A-F of Figure 1 illustrated that the integration of carrier phase observation and pseudorange led to an IP with a higher magnitude precision, and the performance of the methods improved with the precision of RP. The noise level of the smoothed pseudorange obtained from Hatch's filtering was decreasing by expanding the smooth window. The noise seemed unchanged as we were putting more epochs' data into the filter, as the filter could only provide a limited precision. The comparison results among various LS methods showed that there was little difference among the DELS, TELS, and QELS, which was difficult to distinguish in the subgraph A-F of Figure 1; a detailed comparison is presented in the following section. The comparison results obtained from different types of methods showed that the LS strategy could lead to enhanced performance with respect to Hatch's method, even for k = 100, which was becoming more and more prominent as enlarging the noise level of pseudorange, especially at the situation where the precision of pseudorange is 3 m and Hatch's method obtained the lowest accuracy at the test. Besides, Figure 1 shows that the LS required some epochs for convergence, while it was not required for the Hatch's methods because the Hatch's sacrificed some observation information. For instance, Hatch's method could not provide smoothed pseudorange for the first 20 epochs if the k = 20, for the 40 epochs if the k = 40, etc. Figure 1 reveals the residual error comparison curve of the Hatch's smoother and the types LS method, which might verify that the LS method performed better than Hatch's filter at any precision pseudorange. However, to further evaluate the performance of the two types of methods, we analyzed their performance using the RMS (Root Mean Square) as the criterion, for different epoch numbers participating in the filter and different noise level of RP, expressed as std (standard deviation), presented in Table 1. Due to the insignificant difference among the various LS methods, such as DELS, TELS, and QELS, only DELS was used in this analysis. Additionally, in this analysis, the first several epochs were not included, where the LS method produced relatively sizeable residual error IP, which would significantly expand the RMS of the LS method but not for Hatch's smoother.
IGS organization. The emulator outputted measurements at DOY (day of the year) 300 (27 October), 2018 for ALBH station, which was located in Canada. Based on the reference provided and satelliteto-receiver distance mentioned above, it took the E02 (the GALILEO navigation satellite with pseudorandom noise code of 2) as an example to test and verify the performance of Hatch's smoother, the DELS, and the M/TELS method. The subgraph A-F of Figure 1 shows the residual error of various methods as changing the precision of pseudorange accuracy, where the k represents the smoothing time constant of Hatch's smoother. The residual error presented in subgraph A-F of Figure 1 illustrated that the integration of carrier phase observation and pseudorange led to an IP with a higher magnitude precision, and the performance of the methods improved with the precision of RP. The noise level of the smoothed pseudorange obtained from Hatch's filtering was decreasing by expanding the smooth window. The noise seemed unchanged as we were putting more epochs' data into the filter, as the filter could only provide a limited precision. The comparison results among various LS methods showed that there was little difference among the DELS, TELS, and QELS, which was difficult to distinguish in the subgraph A-F of Figure 1; a detailed comparison is presented in the following section. The The comparison between Hatch's filters and LS methods showed that the performance provided by the two types of filters became more prominent when the noise of RP was enlarged. Overall, the Hatch's filter could reduce the noise of pseudorange by order of magnitude or more, which provided a smoothed pseudorange with tens of centimeters of noise if the precision of RP was 1 m or even several centimeters of noise for the precision of RP of 0.1 m or 0.3 m. Besides, the performance of Hatch's filter was improving with the enlargement of the smoothing windows, such as the performance of smooth windows with k = 100 was often better than k = 80 under various accuracies of RP, which was consistent with the theoretical analysis results provided by Equation (8). Moreover, the filters seemed to tend to stabilize at breakneck speed if it provided a sizeable smooth window or smoothing time constant. Hatch's filter obtained a better performance with a larger smoothing window, but it sacrificed more observation information. However, the RMS of the DELS method was smaller than the Hatch's method at any epoch with various precisions of RP. The RMS of the LS method also seemed to be smaller, with the increasing epoch participating in filtering. Additionally, the detailed comparison displayed that the DELS had a little difference comparing with Hatch's smoother at k = 80/100, where the noise level of pseudorange was 0.1 m or 0.3 m, confirming that the LS method obtained indistinguishable advantage comparing to Hatch's smoother at a situation where pseudorange provided high precision or close to carrier phase precision.
Furthermore, it was proved that the RMS increased as pseudorange accuracy decreased, and the convergence speed of the LS method reduced as the epoch number participating in the filtering expanded in most cases. Theoretically, the LS method may achieve better performance than Hatch's smoother for integrating pseudorange and carrier phase observation if receivers continuously receive signals from the corresponding satellites for long-times. Moreover, this study also investigated the filtering epoch number required to achieve a specific percentage of pseudorange accuracy to compare the convergence speed of the two methods, which is illustrated in Table 2. The convergence defined by other studies is obtaining residual error less than the predefined threshold at the current epoch and the following twenty measurements [33,34]. The event on this question is that the residuals error, of 20 consecutive residuals error, from the integrated pseudorange, less than the threshold, makes it converged to the value.  Table 2 detailly illustrates that a relatively small window of Hatch's filter required more epochs than a relatively larger window in order to achieve an absolute percentage precision of the RP under any situation, which showed that it could achieve higher filtering accuracy with the increasing filtering window and sacrificing more observation information. The required epochs of Hatch's smoother seemed to increase with decreasing the precision of pseudorange, but it made no difference for the LS method. By decreasing the percentage of filtering accuracy in pseudorange accuracy, the number of epochs required to achieve an absolute precision increased, meaning that the speed of the filters' convergence decreased and would converge to a particular value, and it proved that both filters belonged to the convergence filters. However, the LS method might achieve better performance than Hatch's smoother even if the smoothing window was 100 (k = 100) under various pseudorange accuracies. The LS method could achieve 5% precision of RP, although it required a lot of epochs, which might not be reached by the Hatch's method, such as the precision of RP was 0.1 m and the k = 20/40, etc. Several conclusions were derived from the analysis of the comparison between the methods of integrating the carrier phase and pseudorange by using simulated data. However, those results were affected by some other reasons, and thus the types of Hatch's filters and LS were compared based on the theoretical accuracy in this study, which was achieved under various precision of RP. The theoretical accuracy of the tested integrating approaches mentioned was based on the adding 3 mm noise level to carrier phase observation and various noise levels to the RP (Figure 2). To give a bright contrast between the two methods, it only expressed the solution for 1000 epochs caused by the fast divergence of the LS method.  Figure 2 provides the theoretical accuracy of the Hatch's smoothers and LS methods based on various smoothing windows, which theoretically explained the reasons for the phenomena in Figure  1 and Table 2. Based on Figure 2, it was observed that (i) the theoretical filtering precision of Hatch's smoother improved when the smoothing window for the noise level of RP was enlarged, (ii) the noise level of IP at any epoch was higher than the previous epochs, which was consistent with the Equations 7-11, (iii) the accuracy of the initial smoothed pseudorange was equal to the LS method, which proved that the approach to obtain the initial value for Hatch's method was a type of LS method, (iv) the theoretical accuracy of any Hatch's method was less than the LS method, even if k = 100 or more abundant, (v) the precision of smoothed pseudorange was close to LS method as we enlarged the smoothing window, which explained the phenomenon displayed by the simulation data and the result obtained from Table 1.

The Comparison of Least Squares Method Based on Various Epochs
From the comparison between the various methods of Hatch and LS methods, the enhanced performance of LS methods against Hatch's filter was revealed. However, the performance of various LS methods had not been analyzed due to the small differences between the methods. In Table 3, the RMS derivation from the analysis of DELS, TELS, and QELS is presented.   Figure 1 and Table 2. Based on Figure 2, it was observed that (i) the theoretical filtering precision of Hatch's smoother improved when the smoothing window for the noise level of RP was enlarged, (ii) the noise level of IP at any epoch was higher than the previous epochs, which was consistent with the Equations (7)-(11), (iii) the accuracy of the initial smoothed pseudorange was equal to the LS method, which proved that the approach to obtain the initial value for Hatch's method was a type of LS method, (iv) the theoretical accuracy of any Hatch's method was less than the LS method, even if k = 100 or more abundant, (v) the precision of smoothed pseudorange was close to LS method as we enlarged the smoothing window, which explained the phenomenon displayed by the simulation data and the result obtained from Table 1.

The Comparison of Least Squares Method Based on Various Epochs
From the comparison between the various methods of Hatch and LS methods, the enhanced performance of LS methods against Hatch's filter was revealed. However, the performance of various LS methods had not been analyzed due to the small differences between the methods. In Table 3, the RMS derivation from the analysis of DELS, TELS, and QELS is presented.  Table 3 gives the RMS error for different LS methods under various precisions of pseudorange, which showed that the accuracy for all the methods generally decreased as the precision of RP reduced. The rate of RMS decrease was gradually reduced with the increase of the number of epochs used in the solution. In some of the examined cases, the solution obtained from the TELS/QELS methods was better than the DELS method. However, the differences between the three methods were generally very marginal, making uncertain any conclusions based on the simulations. Thus, the three methods were further evaluated based on their theoretical accuracy as it was calculated as a ratio according to the error propagation theorem.
Subgraphs of A-F in Figure 3 confirms that there were tiny differences among different types of LS method, especially the difference between TELS and QELS, as changing the precision of pseudorange. Besides, the gap would be further reduced if we increased the number of epochs participating filter at one time, which would be caused by enlarging the self-correlation with increasing epoch number at one-time filtering. There existed 2 non-zero elements for the DELS, 8 for the TELS, and 18 for the QELS at the off-diagonal of the variance-covariance matrix for the LS method if the filter entered the loop entirely. Additionally, it was difficult to obtain better accuracy for the smoothed pseudorange if the precision of pseudorange was low, even if it had been filtering for tens thousands of times, such as the noise level of the RP was more significant than 0.8 m. Thus, the accuracy of the LS methods were decreasing as we decreased the precision of the original observation information. The MELS might get more obvious accuracy improvement for the GNSS devices with high accuracy and providing a long-time observation. Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 20

The Real Measured Data Test Results and Performance Analysis
The proposed method was also evaluated by using real GNSS data. A GNSS receiver was mounted on a vehicle recording the positioning of the vehicle, which was moving on a predefined trajectory on the athletics field of China University of Mining and Technology (Figure 4). The GNSS receivers recorded data on DOY 278 (5 October) of 2019, continuing for appropriately 1000 s, with 1 Hz sampling frequency. The effectiveness of the method was assessed based on the positioning accuracy.

The Real Measured Data Test Results and Performance Analysis
The proposed method was also evaluated by using real GNSS data. A GNSS receiver was mounted on a vehicle recording the positioning of the vehicle, which was moving on a predefined trajectory on the athletics field of China University of Mining and Technology (Figure 4). The GNSS receivers recorded data on DOY 278 (5 October) of 2019, continuing for appropriately 1000 s, with 1 Hz sampling frequency. The effectiveness of the method was assessed based on the positioning accuracy.
The GNSS receiver recorded eight full laps on the athletic track, and the trajectory had a variety of characteristics and could represent most cases of motion. Its scenario included linear and curvilinear motion, which consisted of acceleration, deceleration, and uniform motion (the stationary state is a particular case of uniform motion). The receiver was at a standstill at the beginning and started to move around 260 s with high acceleration. Then, it moved in a circular motion with a constant speed after accelerating to about 4 m/s. Next, it continued to accelerate and decelerate until it completed the movement and returned to the starting point. The variation of the velocity is expressed in Figure 5, where D indicates the 3D (three directions) speed, and E, N, and U represent the East, North, and Up directions, respectively, whose coordinate origin is the average of all sample points' coordinates. The GNSS receiver recorded eight full laps on the athletic track, and the trajectory had a variety of characteristics and could represent most cases of motion. Its scenario included linear and curvilinear motion, which consisted of acceleration, deceleration, and uniform motion (the stationary state is a particular case of uniform motion). The receiver was at a standstill at the beginning and started to move around 260 seconds with high acceleration. Then, it moved in a circular motion with a constant speed after accelerating to about 4 m/s. Next, it continued to accelerate and decelerate until it completed the movement and returned to the starting point. The variation of the velocity is expressed in Figure 5, where D indicates the 3D (three directions) speed, and E, N, and U represent the East, North, and Up directions, respectively, whose coordinate origin is the average of all sample points' coordinates. Before we processed the GNSS data, the software package loaded the orbit and clock difference information and then detected and repaired the cycle slip by the GF (Geometry-Free) combination [35] and MW (Melbourne-Wubbena) combination [36,37]. The cutoff of the elevation angle was set to 10° to decrease the error on the propagation path. This study eliminated the ionospheric delay by IF combination [38], compensated for the static delay of tropospheric delay [39], and ignored dynamics delay. Additionally, it compensated and transformed the related errors into the corresponding satellite-earth distance on the signal propagation direction based on the relevant model. The model estimated the clock error and position of the receiver as parameters based on the LS. The accurate trajectory of this analysis was by analyzing the GNSS data with IE (Inertial Explorer), a postprocessing software developed by the Waypoint product group of the NovAtel, which processed the GNSS datasets by PPP technology with precision products provided by GFZ (Helmholtz-Centre Potsdam -German Research Centre for Geosciences). Based on the high-accuracy positioning result, Table 4 analyzed the true error of various methods and calculated the statistical characteristics (mean, std, and RMS) transformed into E-N-U coordinate system (Table 4).  The GNSS receiver recorded eight full laps on the athletic track, and the trajectory had a variety of characteristics and could represent most cases of motion. Its scenario included linear and curvilinear motion, which consisted of acceleration, deceleration, and uniform motion (the stationary state is a particular case of uniform motion). The receiver was at a standstill at the beginning and started to move around 260 seconds with high acceleration. Then, it moved in a circular motion with a constant speed after accelerating to about 4 m/s. Next, it continued to accelerate and decelerate until it completed the movement and returned to the starting point. The variation of the velocity is expressed in Figure 5, where D indicates the 3D (three directions) speed, and E, N, and U represent the East, North, and Up directions, respectively, whose coordinate origin is the average of all sample points' coordinates. Before we processed the GNSS data, the software package loaded the orbit and clock difference information and then detected and repaired the cycle slip by the GF (Geometry-Free) combination [35] and MW (Melbourne-Wubbena) combination [36,37]. The cutoff of the elevation angle was set to 10° to decrease the error on the propagation path. This study eliminated the ionospheric delay by IF combination [38], compensated for the static delay of tropospheric delay [39], and ignored dynamics delay. Additionally, it compensated and transformed the related errors into the corresponding satellite-earth distance on the signal propagation direction based on the relevant model. The model estimated the clock error and position of the receiver as parameters based on the LS. The accurate trajectory of this analysis was by analyzing the GNSS data with IE (Inertial Explorer), a postprocessing software developed by the Waypoint product group of the NovAtel, which processed the GNSS datasets by PPP technology with precision products provided by GFZ (Helmholtz-Centre Potsdam -German Research Centre for Geosciences). Based on the high-accuracy positioning result, Table 4 analyzed the true error of various methods and calculated the statistical characteristics (mean, std, and RMS) transformed into E-N-U coordinate system (Table 4). Before we processed the GNSS data, the software package loaded the orbit and clock difference information and then detected and repaired the cycle slip by the GF (Geometry-Free) combination [35] and MW (Melbourne-Wubbena) combination [36,37]. The cutoff of the elevation angle was set to 10 • to decrease the error on the propagation path. This study eliminated the ionospheric delay by IF combination [38], compensated for the static delay of tropospheric delay [39], and ignored dynamics delay. Additionally, it compensated and transformed the related errors into the corresponding satellite-earth distance on the signal propagation direction based on the relevant model. The model estimated the clock error and position of the receiver as parameters based on the LS. The accurate trajectory of this analysis was by analyzing the GNSS data with IE (Inertial Explorer), a post-processing software developed by the Waypoint product group of the NovAtel, which processed the GNSS datasets by PPP technology with precision products provided by GFZ (Helmholtz-Centre Potsdam -German Research Centre for Geosciences). Based on the high-accuracy positioning result, Table 4 analyzed the true error of various methods and calculated the statistical characteristics (mean, std, and RMS) transformed into E-N-U coordinate system ( Table 4).
The comparison results among the Hatch's methods from Table 4 showed that the mean value of E, U, and 3D directions increased, while the N direction decreased, and the std of all directions decreased as increasing the value of k, proving that enlarging the smoothing window would mainly reduce the positioning noise and then the mean. Moreover, it was shown that both Hatch's filter and LS method could improve the pseudorange accuracy, thereby the positioning accuracy, with providing a smaller RMS value than the RP solution, but the Hatch's filter and LS method enhanced the positioning performance by different methods. As we compared the positioning result obtained from the RP and Hatch's filter, it showed that the mean of E, U and 3D directions of Hatch's filter was more extensive than RP in most cases, while the mean of N direction and the std of positioning result provided by Hatch's methods at all directions was smaller than RP. The comparison results proved that the smoothed pseudorange improvement of the positioning accuracy firstly relied on reducing the noise and then the mean value at N direction. As we compared the positioning result obtained from the RP and LS filters, it showed that the mean value at E, U and 3D directions, and the std at all directions obtained from the types of LS method were smaller than RP, which proved that the LS method firstly enhanced the performance by reducing the std and then the mean value at E, U and 3D directions. Additionally, Table 5 summarizes the percentage accuracy improvement of various methods compared to solution of RP. Table 5 shows that the types of LS method improved about 28% at E, 6% at N, and 12-13% at U direction, respectively, while the Hatch's methods improved about 9-16% at E, 9-14% at N, and 7-9% at U direction, respectively. However, the types of LS method and Hatch's methods improved about 12% and 9-10% at the 3D direction, respectively, which illustrated the accumulating accuracy improvement at the E and U directions comparing LS method with Hatch's methods was more than the decreased at the N direction, proving that the types of LS method could obtain a better performance than Hatch's method. As we compared the different types of LS methods, the TELS method could provide a little better performance than DELS while the QELS could provide a little higher positioning accuracy than TELS, proving that the LS methods based on multi-epochs only provided a small difference with several millimeters corresponding to several percentage points, which is caused by the self-correlation of the constructed observation matrix.

Conclusions
To improve the precision of carrier phase smoothing pseudorange and expand the application of GNSS, we enhanced the DELS to TELS, as well as its extension method named MELS, for integrating the pseudorange and carrier phase observation, which relied on three or more epochs. To compare the performance of the developed methods against the Hatch's smoother, the methods were analyzed by using simulated data of various precisions of RP. The numerical results from satisficing the residual error showed that both two types of filtering methods belonged to convergence filter. The test result from simulation shown that the convergence accuracy and speed were analyzed and compared theoretically, which showed that the types of LS could have a better performance than Hatch's methods. The experiment results from measured GNSS data shown that the types of LS provided better performance than Hatch's method at E and U directions. Moreover, the types of LS method and Hatch's methods improved about 12% and 9-10% at the 3D direction, respectively, which illustrated the accumulating improvement at the E and U directions comparing the LS method with Hatch's methods was more than the decreased at the N direction, proving that the types of LS method could obtain a better performance than Hatch's method. Additionally, the DELS, TELS, and QELS were compared in this study, and it was revealed that the TELS method could provide a little better performance than DELS while the QELS could provide a little higher positioning accuracy than TELS at 3D direction, proving that the LS method based on multi-epochs achieved only a millimeter accuracy difference compared to that of the LS method based on two epochs.