A Novel Dual-Domain Filtering Method to Improve GNSS Performance Based on a Dynamic Model Constructed by TDCP

To overcome the strikes from the tremendous noise of pseudorange and the integer ambiguity of carrier phase observation, we developed a novel dual-domain filtering method that integrates the traditional observation- and an innovative position- domain filtering method. The novel method employed the conventional Hatch’s smoother to improve the precision of pseudorange at the observation-domain filtering part. At the position- domain filtering part, the proposed method constructs an accurate dynamical model for KF (Kalman filter) or ARKF (adaptive-robust KF) that relies on the TDCP (time differential carrier phase) velocity measurement technology. To assess the effectiveness of the novel method, we compared it with other methods under various system integration by measured data from IGS (International GNSS Service) MGEX (Multi-GNSS Experiment) and dynamic data collected at CUMT (China University of Mining and Technology). The static numerical experiments obtained from IGS datasets show that the ARKF improves the precision of 60/84 and 55/84 situations’ directions compared with the LS method based on OP (original pseudorange) and SP (smoothed pseudorange), respectively, and Hatch’s smother improves the precision of all situations’ directions. The experiment results prove that both types of observation- and position- domain filtering methods can improve positioning precision. Also, the novel dual-domain filtering improved the precision of 80/84 situations’ directions compared with the LS method. Moreover, the dynamic experiment shows that the novel method can improve all of the situations’ directions except the N (North) direction provided by KF base on SP under GC and U (Upon) direction provided by LS base on SP under GR. The average of improved accuracy is 0.212m, 0.346m, and 0.588m at the E (East), N, U direction and its corresponding average percentage of improved accuracy is 30.2%, 60.2%, and 26.7%, respectively, which proves the novel dual-domain filtering method provided the best performance among the tested algorithms.


I. INTRODUCTION
After GPS (Global Positioning System) and GLONASS (GLObalnaja NAwigazionnaja Sputnikowaja Sistema) providing global services, the BDS (BeiDou Navigation Satellite System) II has achieved comprehensive network in the Asia-Pacific region with launching 14 satellites by 2012 and BDS III planned to complete the global networking by 2020 [1]. At the same time, the GALILEO (GALILEO Satellite Navigation System) designed for European civilians has The associate editor coordinating the review of this manuscript and approving it for publication was Larbi Boubchir . launched a series of satellites. The GNSS has an enormous influence and plays an irreplaceable role in human's daily life, such as military, national defense, transportation, agriculture, forestry, and other industries.
However, the strikes from the tremendous noise of pseudorange and the integer ambiguity of carrier phase observation significantly limit the positioning accuracy and application of GNSS. Thus, generating new measurements that provide high-precision and non-existent ambiguity is a crucial issue to improve positioning accuracy. In general, it can divide the existing filters into observation-and position-domain method [2], where the Hatch's method is the first proposed and widely used in the observation-domain [3]. Pieces of literature are devoted to enhancing filters by optimizing various aspects, such as ionospheric variation [4]- [8], multipath effects [9], [10], the smoothing window [11], the correlation of current and previous [12]. Also, some scholars soothed the pseudorange algorithm based on the parameter estimation theory [2], [13]- [15] or doppler observation [16]- [19]. Another type of filtering method is the Kalman filter (KF) that effectively utilizes historical measurement information and has been studied by experts and scholars [20]. Yang et al. proposed an IGG (Institute of Geodesy and Geophysics) III function solution to construct the equivalent weight for the robust KF [21], and Guo restrains the influence of gross error based on a new adaptive factor that treats the significance of observation as a whole [22]. Then, Bahadur et al. obtained a reliable solution relying on the model provided by Guo and developed a PPP (precise point position) software package called PPPH [23]. In general, the position-domain filters obtained better performance than the observation-domain filtering at the same measurement condition [12], [24], but few studies applied both of the models to position and navigation. To improve the positioning accuracy of GNSS and extend its application, we proposed a novel dual-domain filtering method that integrates the traditional observationand an innovative position-domain filtering method in this paper.
The dual-domain filtering method employed Hatch's method to decrees the noise of pseudorange at the observation-domain filtering part [3]. At the position-domain filtering part, the proposed method constructed an accurate dynamical model for KF [20] and ARKF (adaptive-robust KF) [22] that bases on the TDCP (time differential carrier phase) [25].
In the following part, we introduced the observation, stochastic, and dynamic models firstly, and then deduced two types of filtering methods employed by the novel method. Next, this manuscript tested the static datasets from IGS (International GNSS Service) organization and dynamical receiver at CUMT (China University of Mining and Technology) under various system integration, such as G (GPS), GC (GPS/BDS), GR(GPS/GLO), and GRC (GPS/GLO/BDS). The last part made a summary of the contribution and limitation of the proposed method.

II. MATHEMATICAL MODEL OF MULTI-GNSS
A. OBSERVATION MODEL Models should consider the features of different systems for integrating multi-GNSS, including coordinate reference frame, time starting point, signal frequencies. Thus, we applied the ITRF (International Terrestrial Reference Frame) coordinate frame and GPST (GPS time) system obtained from IGS products in this paper. The organization provides post-processing, real-time, predicted satellite orbits, and clock error products in the same coordinate and time system after 1994.
Where the superscript s and j represents GNSS system label and satellite index that includes GPS, GLONASS, GALILEO, and BDS; The superscripts r and i represent the receiver and the signal frequency, respectively; The P and ϕ represent pseudorange and carrier phase observations, respectively; λ represents the signal wavelength, and c is the speed of light in vacuum; The τ s , τ r express the clock offset of satellite and receiver; τ s is the time bias of the corresponding system to GPS time; The ρ represents geometric distance as a function of the receiver and satellite coordinates; The I and T represent ionospheric delay and tropospheric delay, respectively; The b and d are hardware delay on pseudorange caused by satellite and receiver, respectively; The B and D are hardware delay on carrier phase observation caused by satellite and receiver, respectively; The δρ represents the multipath effect delay and ε is the residual error of corresponding observations. The products provided by GFZ (Helmholtz-Centre Potsdam-German Research Centre Geosciences) eliminated the satellite orbit error and the clock difference error in the Equations (1) and (2). More, the SPP (standard point positioning) model of multi-GNSS neglected the hardware delay error and corrected the troposphere delay by the Saastamoinen model [26] based on the meteorological date from the global model. Also, the model used in this paper compensated other modellable errors according to the corresponding model. Thus, the multi-GNSS ionospheric-free linear combination can take a form as follows: where the IF is the ionospheric-free linear combination of pseudorange observation, while the ε IF delegates residual error. The f is the signal frequency.

B. STOCHASTIC MODEL
An accurate positioning model includes not only a functional relationship but also a precise stochastic model that correctly reflects the accuracy of observation accurately. We take a unit weight for GPS and GLONASS, and a lower for BDS and GALILEO because of their performance [23]. Thus, the weight ratio of all pseudorange measurements takes form as follows: σ g : σ r : σ e : σ c = 1 : 1 : 2 : 2 It believes a strong relationship between the signal precision and transmission path, and thus determined the VOLUME 8, 2020 stochastic model for each satellite by trigonometric function model that takes modus as follows: Then, the error propagation theorem is applied to obtain the stochastic model for the ionospheric-free linear combination observation, which expresses as follows:

C. DYNAMIC MODEL
There are three methods to obtain velocimetry in the GNSS measurement model, which are doppler, position difference, and TDCP technology. This manuscript constructs a dynamic system model by TDCP velocimetry technology because the technology provides velocity with millimeter/s-level accuracy. Under the model that it only takes the velocity and clock drift in consideration when estimating parameter in order to increase the degree of freedom, the TDCP technology takes a form as follows [25]: where the t m (m = 1, 2) is the signal receiving time. The p s is the satellite position obtained from satellite ephemeris. The p r is the receiver position that can get from Equation (1) based on SPP technology. The e(t m ) is directional cosine coefficient decided by the receiver and satellite coordinates . The ∇ t r is the receiver clock drift, which is also called clock velocity.
If it sets ζ = e(t 2 )· p s (t 2 )− e(t 1 )· p s (t 1 )+ p r (t 1 )[ e(t 1 )− e(t 2 )], it can rewrite Equation (7) as follows: where V r indicates the receiver velocity. According to model (8), the one-step prediction of the carrier's position and its variance-covariance as follows: where t represents the sampling interval; δ is the std (standard deviation) of the corresponding estimator.

III. FILTERING METHOD
The observation-and position-domain filters are two categories of filtering methods in GNSS technology, and this section introduces the novel dual-domain filtering method applied in this paper.

A. POSITION-DOMAIN FILTERING METHOD
A given stacking system expresses its state-space model as follows: whereX k andX k−1 represent the state vector of the system at the epoch k and k − 1, respectively; Z k is the measurement signal at the current time; F k/k−1 is the state transition matrix; k/k−1 is the noise allocation matrix; H k is the measurement matrix; W k−1 is the system noise vector and V k the measurement noise vector. The unrelated W k−1 and V k are both gauss white noise vector sequences and the corresponding covariance is Q k and R k respectively.
Studies proposed various ARKF methods to get a balance between the one-step prediction error and measurement error, one of which relies on the concept of equivalent weight with a hazard function as following [21], [27]: (12) where theR k is the robust equivalent weight matrix of the Z k which is the adaptive estimation of the observation vector, the α k (0 < α k ≤ 1) is the adaptive factor. TheX k,k−1 is the one-step prediction of the state vector realized by the solution from Equation (9). TheX k and P k is the variance matrix of the updated state vector. This model set 10m for the precision of SPP at a full GNSS environment, and a recursive formula of the filtering method is illustrated as Equations (13)- (17).
Depend on the dynamic model obtaining from TDCP technology, one-step prediction of the state vector,X k,k−1 , takes a form as:X Based on error propagation theory, the covariance matrix, P k,k−1 , of one-step prediction of state vector using a dynamic model: According to the free extremum principle and Lagrange multiplier theory, the gain matrix, K k , is expressed as following that reflects the relationship between the predicted vector and measurements: Thus, the state estimation,X k , and its variance-covariance matrix, P k , at the filtering epoch with a form as following: From the above formulas, the difference between the KF and ARKF is the value of α k . If α k = 1, it would be the conventional extended KF. If 0 < α k < 1, it would be a robust solution with reducing the precision of one-step prediction. Thus, the adaptive factor should appropriately reflect the precision and error relationship between measurements and one-step prediction. One of the methods that treat the influence of observation as a whole to obtain the adaptive element describes as follows [22]: where theṽ k,k−1 is the residual observation error that gets from the predicted state vector with a valueṽ k, The c 0 and c 1 are two constants with experienced values, c 0 = 1.5 ∼ 3.0 and c 1 = 3.0 ∼ 8.0. We give the parameters that c 0 = 2.5 and c 1 = 6.5 in this paper. The n k is the number of observations. The α k is the adaptive factor with a value between 0 and 1 that controls the whole weight ratio of the new measurement information.

B. OBSERVATION-DOMAIN FILTERING METHOD
The observation-domain filtering method is improving the precision of the original measurement based on the smoothing algorithm before parameter estimation. This paper used an enhanced carrier smoothing pseudorange method that considered the ionospheric variation to improve the precision of the SP (smoothed pseudorange). For convenience, we neglected the other error except for ionospheric delay when deducing formula. The Hatch's smoothing method generates a new observation with high precision and nonexistent integer ambiguity based on the advantages of carrier phase and pseudorange that is realized by the following Equation [3]: where S is the SP by carrier phase observation, and M is the smoothing coefficient as we are smoothing pseudorange by carrier phase observation.
Equation (20) clearly shows that the receiver should continuously lock the carrier phase observation. If the receiver lost lock or cycle slip occurs, the integer ambiguity would make a change, and then the smoother should be reset. Thus, the GF (geometry-free) combination and M-W combination [28]- [30] model was employed to detect and repair cycle slip to ensure smoother working continuously. There is another problem that Equation (20) ignores the variety of ionospheric delays and believes the ionospheric delay is a constant, which is always unsuitable for real environments, especially at conditions with disturbing ionospheric activity or the smoother working for a long time. However, the variety of ionospheric delay makes no sense if we use the IF combination for positioning because it can eliminate the value of the error accumulated and absorbed.
Maybe another critical problem of Hatch's smoother is how to determine the initial SP. If there exists a significant deviation on S s,j i,r (t 0 ), the smoother may run for a long time to weaken the systemic error. Thus, it should provide as accurate as possible, and one of the optional methods is employing several epochs to obtain the initial SP, making a mean of the value and providing a model as follows: where the t 0 is the first epoch as long as the GNSS signal is locked. The k is the epochs used to calculate the initial value of SP. The n is a variable.

IV. EXPERIMENTAL ANALYSIS
This section makes a description of the datasets and introduces the datasets processing method, and used the static datasets from IGS MGEX organization and the dynamic data collecting at CUMT to compare the novel method with others under various system integration.

A. DESCRIPTION SOLUTIONS
To thoroughly test the performance of the novel dual-domain filtering method, we compared the proposed method with others by static and dynamic datasets, which bases on position-domain filters with/without the observation-domain filtering under various system integration. The positiondomain filtering methods include the LS (least square), KF, and ARKF, and the system integrations contain a single-system of GPS, dual-system of GC and GR combinations, a triple-system of GRC combination. The GNSS data processing strategy is as follows: (i) the program loads the files of the satellite orbit, clock difference (provided by precision products or broadcast ephemeris), and attitude information of antenna. (ii) the package detected and repaired cycle slips based on the MW and GF linear combination [28]- [30], and checked the continuity of the clock [31]. (iii) 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, KF, and ARKF with a dynamic constructed by TDCP.

B. VERIFICATION RESULTS OF STATIC DATA
Considering the coverage area of BDS, we evaluated the novel double-domain filtering method relying on measured VOLUME 8, 2020 datasets of Asia-Pacific provided by IGS MGEX. Since there may be differences and influences caused by constructing stations and the type of the receiver or antenna, we also tested the algorithm by a dataset from Hong Kong CORS (Continuously Operating Reference Stations). All of the datasets are all collected with 1Hz sampling frequency from 10:00 to 11:00 am on doy (day of the year) 200 (July 19), 2019. Figure 1 shows the distribution of all stations. The distribution of stations shows the static datasets have good representativeness and diversities since the datasets are collected all over the Asia-pacific region. From the distribution of the north to south direction, there are 4 stations in the northern hemisphere, 3 stations in the southern hemisphere with 1 station nearby equator. From the distribution of the east to west direction, 2 stations are standing the west and 5 stations in the east with 90 • E as the dividing line. Moreover, the datasets from inland and littoral are both selected in the datasets. Then, we processed the static datasets by simulating dynamic data, which says that its processing strategy used is the same as the dynamic data. To thoroughly analyze the effectiveness of the proposed approach, we estimate the solutions by methods, such as the least square, KF, and ARKF with OP and SP. This paper counted the RMS (Root Mean Square) of the tested methods under various system integrations for accuracy analysis is a list as the following formula, and Figures 2-4 show the statistical result: where V is the real error based on an actual reference value provided by MGEX organization and IE (Inertial Explorer) software developed by Novatel. n_epoch is all of the epochs at the project.
The static datasets experiment result shows that the accuracy in the E and N direction is higher than U direction, and the multi-GNSS has a better performance than a single-or dual-system in most cases while the double is better than the single, and the conclusions are consistent with the typical   conclusion. Moreover, most of the positioning accuracy of GR is higher than GC, which is proved by most of the stations. The CUT0 is the most special station because there may be a gross error on the GLONASS observation that causes the KF to diverge when estimating parameters. Therefore, it leads to significant deviations on the solution for the ARKF based on the OP (original pseudorange) and KF based on SP under the GR system combinations. Next, we investigate the three-direction accuracy of the testing stations based on various algorithms when comparing the novel with other algorithms. The GNSS performance compared result between LS and ARKF with a dynamic model constructed by TDCP shows the novel method improves the precision of 60/84 situation's directions based on OP and 55/84 based on SP, which proves that the position-domain enhanced by the TDCP was constructive for GNSS performance. To investigate the performance of the novel dual-domain filtering, we made a comparison between the ARKF method enhanced by TDCP based on the SP and LS method based on the OP, the performance of 80/84 stations' direction has improved. The test result mentioned above proves that the novel method has the best performance among the tested algorithms. Then, we investigated the influence of robustification by comparing the KF and ARKF, and the result shows that the ARKF can improve 76/84 situations' directions based on OP and enhance 65/84 situations' directions based on the SP. However, the difference between the two is about millimeters to centimeters if there is no gross error at the observations, and the result shows that there have a very close performance between the two in a normal situation. However, the novel dual-domain filtering method improves the performance of GNSS significantly compared to other tested approaches, despite there being gross errors on the observations, which has been proved by the test result of CUT0 station. The method of the KF based on SP and the ARKF method based on OP both are failures to obtain an acceptable position solution under system integration of GPS/GLO, while the novel dual-domain filtering can obtain a typical result even if there is a gross error on observation.

C. VERIFICATION RESULTS OF DYNAMIC DATA
Furthermore, to provide a more objective evaluation result of the novel's performance, the field data collecting at the campus of CUMT in Xuzhou, Jiangsu province, China, is used to test the algorithm. The receiver used to obtain dynamic GNSS data is the Trimble R10 that receives dual-or triplefrequencies signals from GNSS satellites such as GPS, BDS, GLONASS, Galileo, and others. The dynamic GNSS data was acquainted at a vehicle on 17 June 2017 with a 1 Hz sampling frequency, continuing for about 930s (it removes some unsolvable epoch at the beginning). Figure 5 shows the trajectory map of the receiver. Similar to the static data, we process the GNSS datasets based on various filtering methods and system integration. Then, the Tables 1-4 shows the RMS of three directions based on corresponding filters and system integration.
From the statistics information of dynamic positioning accuracy, we obtained a similar conclusion to static datasets that is the accuracy in the E and N directions is higher than the U direction, and the result obtained from multi-GNSS has a better performance than a single-or dual-system in most cases. The positioning result of GC combination at E and U direction is better than GR, while the GR is better than GC at N direction, providing a different result from the static datasets. Then, the result compared novel dual-domain filtering method with others shows that the approach improves positioning accuracy of all directions under the strategies   with/without observation-domain filtering. However, there is an exception as comparing the novel with the LS method based on the SP that is the U direction under GPS/GLO system combination. It has a quite close result between the two, whose difference is only 0.023m. Also, the comparison between the observation-domain filtering method illustrated that the performance of SP is better than the OP under any situation. The test result between KF and ARKF with a dynamic model constructed by TDCP shows the two have a very close performance, and the ARKF may be slightly better than KF. The KF has a negligible better performance than ARKF with only 0.004m difference at the N direction based on SP under GPS/GLO system combination. Moreover, Tables 5 and 6 show the comparison of the novel method and other approaches.  As illustrated in Table 5, the novel method improves the position accuracy of 96.7% situation's direction comparing with other methods under various system combinations. The accuracy and percentage improved most is U and N direction based on LS with OP under a single-system of GPS, and the corresponding value is 1.635m and 126.6%. Two situations reducing the positioning accuracy are the N direction based on KF with the observation-domain filtering under GC system combination and the U direction based on LS with the observation-domain filtering under the GR system combination. The reduced accuracy of the two is about 0.004m providing a reduced percentage of 0.6%, and 0.023m providing 0.7%, respectively. Therefore, the experiment result clearly illustrates that the novel method improved position accuracy and provided the best performance among various methods. However, tables 5 and 6 shows that the ARKF method just improves a little or decrease the positioning accuracy compared with KF based on the SP/OP, which indirectly describes that the AR strategy can improve the robustness of the method, necessarily the positioning accuracy. However, there exists some positioning solution from KF that is lower than LS that may because the novel method did not take the correlation for the SP and the TDCP measurements.

V. EXPERIMENTAL ANALYSIS
To improve positioning accuracy and extend the application of GNSS, we proposed a novel dual-domain filtering method that integrates the traditional observation-and an innovative position-domain filtering method in this paper. The proposed method built an accurate dynamic model by the TDCP technology to apply the KF/ARKF method for the SPP technology and decreased the noise of pseudorange by Hatch's smoother filter.
To test the performance of the proposed method, we compare the novel with other filtering methods by the static datasets from IGS organization and the dynamic data collected at CUMT under various system integration, such as a single system of GPS, dual-system of GC and GR, and triple-system of GRC. The test result from IGS MGEX static datasets shows that that the ARKF improves the precision of 60/84 and 55/84 stations' direction compared with the LS method relying on smoothed and OP, and the SP can improve the accuracy of all directions. Also, the novel dual-domain filtering can enhance 80/84 stations' direction compared with LS relying on the OP. Then, the dynamic real GNSS deserts show the novel method can improve positioning accuracy in most cases with an average improved accuracy 0.212m, 0.346m, and 0.5875m at the E (East), N, U directions, whose enhanced rate is 30.2%, 60.2%, and 26.7%. The result above shows that we can get a better GNSS performance by observation-or position-domain filtering method. If the strategy used both the domain filtering method, it could get the highest positioning accuracy among the tested algorithms in most cases. Thus, the dual-domain filtering method with a dynamic model constructed by TDCP had the best GNSS performance among the tested algorithm.
However, the novel method has drawbacks to stick itself, and one of which is that the TDCP technology is only available at the situation where four or more satellites are continuously visible, and the other challenge is that the cycle slip should be detected and repaired. Also, the proposed method ignores the correlation between the SP and the TDCP measurements. In the future, we will introduce doppler observations and construct a rigorous stochastic model to improve positioning accuracy.