Kinematic Precise Point Positioning Using Multi-Constellation Global Navigation Satellite System ( GNSS ) Observations

Multi-constellation global navigation satellite systems (GNSSs) are expected to enhance the capability of precise point positioning (PPP) by improving the positioning accuracy and reducing the convergence time because more satellites will be available. This paper discusses the performance of multi-constellation kinematic PPP based on a multi-constellation kinematic PPP model, Kalman filter and stochastic models. The experimental dataset was collected from the receivers on a vehicle and processed using self-developed software. A comparison of the multi-constellation kinematic PPP and real-time kinematic (RTK) results revealed that the availability, positioning accuracy and convergence performance of the multi-constellation kinematic PPP were all better than those of both global positioning system (GPS)-based PPP and dual-constellation PPP. Multi-constellation kinematic PPP can provide a positioning service with centimetre-level accuracy for dynamic users.


Introduction
Since 1994, the International Global Positioning System (GPS) Service (IGS) organization has provided precise GPS satellite orbit and clock products, enabling the development of a novel positioning methodology known as precise point positioning (PPP) [1,2].Based on the processing of un-differenced pseudo-range and carrier phase observations from a single GPS receiver, positioning solutions with accuracies in the centimetre to decimetre range can be attained globally.Precise point positioning (PPP) is one of the most popular techniques for carrier phase-based precise positioning.PPP sometimes makes use of ionosphere-free linear combination to decrease the effect of ionospheric delay.However, it is not based on integer coefficients, and currently, the state information does not preserve the integer nature of ambiguities.Consequently, PPP cannot adequately resolve ambiguities and access the full range of global navigation satellite system (GNSS) carrier-phase accuracies [3,4].Moreover, long observation times are required for convergence [5].Many researchers have attempted to improve the performance of PPP by improving the precision of the satellite orbit and clock products [2,6] and speeding up the ambiguity-resolution process [7][8][9].Satellite positioning availability and integrity can be significantly improved by using multiple GNSSs so that more satellites will be available [10][11][12].Li et al. tested the accuracy of multiple-constellation PPP and discussed the main challenges associated with this process [13,14].The use of multiple GNSSs is expected to enhance the capability of PPP by improving the positioning accuracy and reducing the convergence time.
Since the International Globalnaya Navigatsionnaya Sputnikovaya Sistema (GLONASS) Experiment (IGEX-98) and the follow-on GLONASS Service Pilot Project (IGLOS), the precise GLONASS orbit and clock data have become available.A combined GPS and GLONASS PPP program was developed by Cai and Gao in 2007 [15].Research on GPS-and GLONASS-based PPP has been conducted continuously since then, and the results have demonstrated that the positioning precision and convergence speed were improved by the dual-constellation signal [16][17][18].The BeiDou Navigation Satellite System (BDS) is a global satellite navigation system that was independently developed, deployed, and operated by China and remains in operation today [19].The BDS consists of two separate satellite constellations: a limited test system that has been operating since 2000 and a full-scale global navigation system that is currently under construction.Table 1 lists the satellites in the BDS constellation at or before December 2012.In total, 23 satellites have been involved, 3 of which are no longer operational.Twenty BDS satellites are currently in operation: 6 in geostationary orbits (GEOs), 8 in 55-degree inclined geosynchronous orbits (IGSOs) and 6 in medium Earth orbits (MEOs).The full constellation is planned to eventually comprise 35 satellites.According to its overall planning schedule, the BDS will have global coverage by 2020 [17].

Common Name
Int. Satellite ID Pseudo Random Noise (PRN) Notes 2012-008A C05 GEO 58.75 Researchers [20][21][22] have also developed a model that combines GPS-and BDS-based PPP.The test results revealed that the combined GPS-and BDS-based PPP can decrease the convergence time and improve the positioning precision.Because of improvements in the precision of the BDS and Galileo satellite orbit and clock products, quad-constellation (GPS, BDS, GLONASS and Galileo) PPP has become possible [23].Tegedor et al. and Cai et al. [24,25] improved the performance of quad-constellation PPP.However, most of the research mentioned above addressed PPP for a statistic object.In this work, the performance of PPP for a kinematic user is addressed.
In this study, we assessed the performance of multi-constellation kinematic PPP in terms of the positioning accuracy and convergence time using the measurements collected from receivers on a vehicle.The multi-constellation kinematic PPP model, Kalman filter and stochastic models applied here are introduced in the second section.Section three describes the multi-constellation kinematic PPP data-processing strategy.The performance of multi-constellation kinematic PPP is elucidated by applying it to real data in the fourth section, which is followed by the conclusions.

Multi-Constellation PPP Model
In this study, we developed a PPP software and conducted several tests of the precision and stability of the data processing achieved using this PPP software.The results showed that the accuracy of statistical positioning could reach the millimetre level in the horizontal direction and the centimetre level in the vertical direction.The root mean square (RMS) value of the tropospheric delay was between 0.01 m and 0.02 m.The convergence time was primarily distributed from 10 min to 40 min.These tests will be described in a future research paper.Compared with other open-source PPP software, this software could batch process large amounts of data, and it used precision products provided by the IGS, Multi-GNSS Experiment (MGEX) and International GNSS Monitoring & Assessment System (iGMAS).The structure design of our software was simpler and had good stability.In this section, the observation equation, Kalman filter model and stochastic model used in our software will be discussed in detail.

Multi-Constellation PPP
Ionosphere-free combination observations are normally used in PPP to remove the first-order ionospheric delay.Their code and carrier phase observation can be expressed as: where f i (i = 1, 2) are carrier-phase frequencies in Hertz; ρ i is the code observation at the ith frequency in metres; φ i is the carrier phase observation in cycles; dt is the receiver clock offset in seconds; dT is the satellite clock offset in seconds; c is the speed of light in metres per second; T trop is the tropospheric delay in metres; N IF are the parameters of the float ambiguity after being redefined; R is the geometric range in metres; ε ρ and ε φ are the code and carrier phase observation noise, respectively, including the multipath in metres; and λ IF is the wavelength after being redefined.
The linearization formula of Equation ( 1) is: In Equation (2), r is the geometric range computed using the linearization point.MF is the mapping function; T ztd is the zenith tropospheric wet delay; ρ IF are the code observations; and Φ IF are carrier observations.x i , y i , z i are the position of satellites, and x p , y p , z p are the coordinates of the stations.
Then, the least squares residuals can be written as: where the parameters dX p , dt, T ztd and N IF are the receiver position, receiver clock offset, zenith tropospheric wet delay and ambiguity of the ionosphere-free combination, respectively; B, A, M, andλ are the corresponding coefficient matrices, respectively.V ρIF and V φIF are the observation residuals.For multi-constellation kinematic PPP, the code hardware delay biases should be absorbed into the receiver clock and system time difference, and the carrier-phase biases related to the frequency should be absorbed into the ambiguity.In addition, the system time difference and receiver clock offset were linearly correlated.Generally, instead of estimating a receiver clock parameter for each satellite system observation, we always introduce a system time difference parameter to reflect the difference between the different system times.The multi-constellation kinematic PPP observation equations can be expressed as follows: where V ρIF and V ΦIF are the GPS observation residuals.V ρIF 1 and V φIF 1 are residuals of other navigation satellites but not GPS satellites.L ρIF and L ρIF 1 are ionosphere-free code observations of GPS and another GNSS, respectively.V φIF and V φIF 1 are the ionosphere-free carrier phases observable.For B 1 , λ other is the corresponding coefficient matrix of other navigation satellites.dt gps is the clock offset of GPS, and dt other is the system time difference parameter.In the tests performed in this study, the ambiguities were treated as constants (assuming no cycle-slips), and other parameters were all epoch dependent, as shown in Equation ( 4).
Combining multiple GNSSs is associated some issues, such as the integration of the time reference system, the coordinate reference systems, and the inter-frequency biases.These issues also represent challenges affecting multiple-constellation PPP.Because the individual GNSSs have different frequencies and signal structures, the code bias values are different in a multiple-GNSS receiver.These biases are the inter-system biases for the code observation.The phase delays are also different, and their differences represent inter-system biases for phase observations.GLONASS satellites emit their signals on individual frequencies, which will also lead to frequency-dependent biases in the receivers.For the GLONASS satellites with different frequency factors, the receiver code biases differ from the carrier phase bias.These differences are usually called inter-frequency biases.Both the inter-system and inter-frequency biases must be considered in a combined analysis of multi-GNSS data.The corresponding parameters should be estimated for all multi-GNSS receivers: one bias for the code measurements of each system and each frequency for GLOANSS.However, the inter-system/inter-frequency biases and obtained satellite clocks are fully correlated.Consequently, when satellite clocks are used, the corresponding biases must also be estimated or corrected for these GNSS receivers.It should be noted that such a receiver internal bias is relevant only when processing the code data.Indeed, when analysing the phase measurements, the corresponding phase ambiguity parameters will absorb the phase delays.These will only be relevant if ambiguities are resolved to their integer values (i.e., mixed ambiguity resolution between different GNSSs, GLONASS ambiguity resolution or un-differenced ambiguity resolution).

Kalman Filter Model
A Kalman filter was used for the multi-constellation kinematic PPP processing in this study.The GNSS dynamic positioning system state equation and observation equation are: where i is the time of t i , and Xi and X i−1 are the m × 1 state vectors at t i and t i−1 , respectively.m is the number of parameters; Φ i,i−1 is an m × m dimensional state transition matrix; and W i is a system noise vector, which is assumed to be drawn from a zero mean multivariate normal distribution with covariance Σ W i .L i is the observation design vector matrix, and e i is the measurement error.P i is a weight matrix.The steps involved in Kalman filtering are as follows [26,27]: (1) Store Xi−1 and ∑ Xi−1 at t i−1 (2) Predicted (a priori) state estimate: (3) Predicted (a priori) estimate covariance: (4) Innovation or measurement residual and covariance: (5) Optimal Kalman gain: (6) Updated (a posteriori) state estimate: (7) Updated (a posteriori) estimate covariance: (8) Let I = i + 1 and then return to the first step until the end of the data.Appropriate stochastic models for the observations and dynamic models for the state vector must be provided in the Kalman filter.The stochastic model is usually defined using an appropriate covariance matrix that describes the statistical properties of the measurements [28].In the following section, the stochastic models used in this study will be introduced in detail.
The kinematic model used in the dynamic model is the constant acceleration: where r(t) is an unknown vector, r(t) is the velocity, and ∆t is the acceleration time length.In our software, the coordinate, receiver clock, tropospheric delay and ambiguity parameters were the unknown parameters.These parameters were estimated using the first-order Gauss-Markov (GM) random process.If no cycle-slips occurred, we treated the ambiguity as a constant.In our software, we only stored the information of the current epoch and the next epoch, and we estimated the parameters of each epoch.The discussion regarding dynamic noise variance of the related parameters is as follows: The discrete first-order GM process is: where x is the state vector, x k = e −∆t/τ , τ is the correlation time, ω is the white noise sequence with zero mean value, and t is the time interval.Generally,β = 1 τ , β is the damping coefficient.If the damping coefficient is too large, the current and next epochs will have greater volatility.However, the two epochs have a strong time correlation.
The state transition matrix is: The dynamic noise variance matrix is: where q is the spectral density or dynamic noise variance matrix.When the correlation time τ is zero, Equation ( 15) represents the pure white noise model.If τ is infinite, the filtering process is a pure random walk.
In general, pure random walks are suitable for simulating the three-dimensional coordinates and ambiguities, where the tropospheric delay and the receiver can be modelled in two forms.
For static PPP without speed and acceleration, the state parameter vector is: The corresponding state transition matrix, Φ k+1,k , and the dynamic noise variance matrix of the three-dimensional position coordinate are: where: q ϕ : latitudinal spectral density; q λ : longitudinal spectral density; q h : elevation spectral density; R m : meridian radius of curvature; R n : radius of curvature on the dome circle; and h: the height of the station.
The dynamic noise variance matrix of the receiver is as follows: If it is a pure random walk: If it is pure white noise: where q dt is the receiver clock density, and β dt is the corresponding damping coefficient.
In general, the tropospheric zenith delay is expressed as a pure random walk, and thus, its dynamic noise variance is: where q trop is the spectral density of the tropospheric zenith wet delay.The ambiguity parameter can be regarded as a constant: Q N = 0.For dynamic PPP, the state parameters should include the speed and acceleration parameters, and the corresponding matrix should change.Thus, the state vector parameters should be: The corresponding state transition matrix is: where: The state noise variance is:

Stochastic Models
Error corrections affect the performance of PPP.Some corrections can be eliminated by using function models, which are discussed in Section 3.However, function models are not sufficient to improve the accuracy and convergence time of PPP.Consequently, an appropriate stochastic model must be applied.Different stochastic models reflect different PPP results.Generally, three common stochastic models exist: the equal-weight stochastic model, the carrier-noise rate-based stochastic model and the elevation angle-based stochastic model.In this study, the last model was applied.
Most of the GNSS observation errors (i.e., the troposphere refraction delay, ionosphere refraction delay, and multipath effect) are related to the elevation angles of satellites.To decrease these errors, stochastic models based on the elevation angles of satellites can be established.Elevation angle-based stochastic models mainly include trigonometric function models and exponential function models [29].In this study, we used the sine function-based elevation angle stochastic model, which is described by Equation (30): where θ is the elevation angle of the satellite, and σ 2 0 is the prior variance of observations.Generally, the multipath and large observation noise exist at low elevation angles.We defined the weight segment to reduce the weight of the observations at lower elevation angles.The corresponding code and carrier phase variance matrices are: where σ 2 P,0 and σ 2 φ,0 are the prior variances of code and carrier phase observations, respectively.α is the elevation angle threshold and is typically set to 30 • .When adopting the code and carrier phase observations simultaneously, the variance-covariance expression is: It should be noted that different GNSSs have different prior observation variances.For the GPS and GLONASS code and carrier phase observations, the precisions are set to 0.3 m and 0.002 m, respectively.Because the BDS satellite orbit and clock have relatively lower accuracies [30,31], their measurements are down-weighted.That is, the phase observation precision is set to 0.004 m, and the code observation precision is set to 0.6 m for the BDS [25].

Multi-Constellation Kinematic PPP Data Processing Strategy
In this study, a self-developed PPP software was used for the data processing.Here, the data preprocessing directly affected the accuracy of the positioning results.The eliminated error and model used in this study are shown in Table 2.The first-order ionospheric delay errors were removed using the ionosphere-free combined observables in our self-developed PPP software.The hydrostatic (dry) tropospheric delay was corrected based on observations using the UNB3m model [32], whereas the non-hydrostatic (wet) part was estimated as a parameter.The Neill mapping functions [35] were used for projection from the slant delays to the zenith delay.The receiver and satellite antenna phase centre offsets (PCO) and phase centre variations (PCV) were collected using the parameters provided by IGS ANTEX.In addition, the MGEX of IGS investigated the new GNSSs [21].In this study, we applied the multi-constellation GNSS precise satellite orbit and clock products from MGEX to mitigate the satellite orbit and clock errors.
The main PPP parameters included the coordinates, the zenith tropospheric delay parameter, multi-constellation system time difference parameters, and the ambiguities.To estimate both the static and kinematic parameters, our software used the Kalman filtering model.The cut-off angle was set at 10 • .To decrease the convergence time and make full use of the satellite observation information, we regarded the satellite as a new satellite when the cycle slip occurred.The cycle slip detection and repair were conducted using the MW combination and geometry-free combination.The software automatically calculated the scaling matrix based on satellite movements to improve the computational efficiency.Simultaneously, the software automatically set the parameters according to the static or dynamic mode.

Experimental Data Description
To assess the availability of multi-constellation kinematic PPP applied to dynamic objects, a multi-constellation kinematic PPP test was conducted in Huainan, China on 18 January 2015.Kinematic GNSS observations were collected over 2 h from three receivers mounted on a vehicle.The speed of the vehicle was approximately 20 km/h.The receivers were Hi-Target V8 (TRM59800.00NONE) devices, which can collect GPS, BDS, and GLONASS observations.The sampling interval was 1 s. Figure 1 shows the vehicle's route, which started from the "start point" and went through each section twice.To investigate the accuracy of the PPP result, we also set a base station with the same type of receiver and antenna at the start point to determine the coordinates of the rover station with cm-level accuracy using the double-difference real-time kinematic (RTK) approach.The surrounding environment was generally good for observations, and the sky visibility was very high.Multi-constellation mixed precise satellite orbit and clock products provided by the Deutsches GeoForschungsZentrum (GFZ) (Germany) were adopted for PPP data processing.It should be noted that no BDS receiver antenna PCO and PCV correction file was available for the antenna mentioned above.Consequently, the BDS PPP positioning results contained a systemic error.

Availability of Multi-Constellation Kinematic PPP
To investigate the benefits of multi-constellation PPP, the availability was analysed in terms of the numbers of visible satellites and the positional dilution of precision (PDOP) value, which is important for quantifying the GNSS positioning accuracy [25].
Figure 2 shows the number of visible satellites and the PDOP values for each processing case.G represents GPS only, G + R represents the GPS/GLONASS combination, G + C represents the GPS/BDS combination, and G + R + C represents the GPS/GLONASS/BDS combination (as mentioned below).As shown in this figure, the GPS-only geometry was weaker than those of the three combined systems.The largest PDOP value of GPS was 20, whereas the PDOP value of the GPS/GLONASS/BDS combination was always below 5.The numbers of satellites and the PDOP values of the four systems varied very stably in the first 0.4 h and from 1.1 to 1.5 h.In contrast, variation occurred very frequently during the rest of the experimental period.The reason for this result was that the vehicle stayed in static mode for the first 0.4 h and the 1.1-1.5 h. Figure 3 also explains this situation, as it shows the coordinate increments at the directions of north, east and up.

Availability of Multi-Constellation Kinematic PPP
To investigate the benefits of multi-constellation PPP, the availability was analysed in terms of the numbers of visible satellites and the positional dilution of precision (PDOP) value, which is important for quantifying the GNSS positioning accuracy [25].
Figure 2 shows the number of visible satellites and the PDOP values for each processing case.G represents GPS only, G + R represents the GPS/GLONASS combination, G + C represents the GPS/BDS combination, and G + R + C represents the GPS/GLONASS/BDS combination (as mentioned below).As shown in this figure, the GPS-only geometry was weaker than those of the three combined systems.The largest PDOP value of GPS was 20, whereas the PDOP value of the GPS/GLONASS/BDS combination was always below 5.The numbers of satellites and the PDOP values of the four systems varied very stably in the first 0.4 h and from 1.1 to 1.5 h.In contrast, variation occurred very frequently during the rest of the experimental period.The reason for this result was that the vehicle stayed in static mode for the first 0.4 h and the 1.1-1.5 h. Figure 3 also explains this situation, as it shows the coordinate increments at the directions of north, east and up.
the numbers of visible satellites and the positional dilution of precision (PDOP) value, which is important for quantifying the GNSS positioning accuracy [25].
Figure 2 shows the number of visible satellites and the PDOP values for each processing case.G represents GPS only, G + R represents the GPS/GLONASS combination, G + C represents the GPS/BDS combination, and G + R + C represents the GPS/GLONASS/BDS combination (as mentioned below).As shown in this figure, the GPS-only geometry was weaker than those of the three combined systems.The largest PDOP value of GPS was 20, whereas the PDOP value of the GPS/GLONASS/BDS combination was always below 5.The numbers of satellites and the PDOP values of the four systems varied very stably in the first 0.4 h and from 1.1 to 1.5 h.In contrast, variation occurred very frequently during the rest of the experimental period.The reason for this result was that the vehicle stayed in static mode for the first 0.4 h and the 1.1-1.5 h. Figure 3 also explains this situation, as it shows the coordinate increments at the directions of north, east and up.

Positioning Accuracy and the RMS of the Kalman Filter
To assess the kinematic positioning accuracy of the multi-constellation PPP, the double difference RTK results were regarded as the true coordinates.These results were computed using the GNSS/INS tightly coupled resolution of Inertial Explorer 8.60 software (IE 8.60), using the multi-constellation signals.The positioning accuracy of this software can reach 1-2 cm. Figure 4a shows the positioning results of multi-constellation kinematic PPP and RTK.When we zoomed in it in Figure 4b, we found that the position of GPS only was very different from that of RTK.There was also little error in the results of the GPS/GLONASS PPP.The positions of the GPS/BDS and GPS/GLONASS/BDS were very similar to those of the RTK.The differences between these RTK coordinates and the PPP results of the GPS only, GPS/GLONASS, GPS/BDS and GPS/GLONASS/BDS were compared to assess the performance of the multi-constellation kinematic PPP in detail.The mean square errors of these differences are shown in Table 3. From Table 3, we can see that the mean square errors of GPS were 0.045 m and 0.085 m in the east and north directions, respectively, which

Positioning Accuracy and the RMS of the Kalman Filter
To assess the kinematic positioning accuracy of the multi-constellation PPP, the double difference RTK results were regarded as the true coordinates.These results were computed using the GNSS/INS tightly coupled resolution of Inertial Explorer 8.60 software (IE 8.60), using the multi-constellation signals.The positioning accuracy of this software can reach 1-2 cm. Figure 4a shows the positioning results of multi-constellation kinematic PPP and RTK.When we zoomed in it in Figure 4b, we found that the position of GPS only was very different from that of RTK.There was also little error in the results of the GPS/GLONASS PPP.The positions of the GPS/BDS and GPS/GLONASS/BDS were very similar to those of the RTK.The differences between these RTK coordinates and the PPP results of the GPS only, GPS/GLONASS, GPS/BDS and GPS/GLONASS/BDS were compared to assess the performance of the multi-constellation kinematic PPP in detail.The mean square errors of these differences are shown in Table 3. From Table 3, we can see that the mean square errors of GPS were 0.045 m and 0.085 m in the east and north directions, respectively, which were much larger than those of GPS/BDS.With the GPS/BDS combination, accuracy improved by 71.11% and 41.18% over GPS only in the north and east directions, respectively.However, the positioning accuracy of GPS/GLONASS was slightly poorer than that of the GPS/BDS.The reason would be that there were more visible GPS/BDS satellites than which of GPS/GLAONSS, and the PDOP of GPS/BDS was better than that of GPS/GLAONSS) (Figure 2).The accuracies of GPS/GLONASS/BDS PPP were better than those of GPS/BDS by 30.77% and 56.00%.

Positioning Accuracy and the RMS of the Kalman Filter
To assess the kinematic positioning accuracy of the multi-constellation PPP, the double difference RTK results were regarded as the true coordinates.These results were computed using the GNSS/INS tightly coupled resolution of Inertial Explorer 8.60 software (IE 8.60), using the multi-constellation signals.The positioning accuracy of this software can reach 1-2 cm. Figure 4a shows the positioning results of multi-constellation kinematic PPP and RTK.When we zoomed in it in Figure 4b, we found that the position of GPS only was very different from that of RTK.There was also little error in the results of the GPS/GLONASS PPP.The positions of the GPS/BDS and GPS/GLONASS/BDS were very similar to those of the RTK.The differences between these RTK coordinates and the PPP results of the GPS only, GPS/GLONASS, GPS/BDS and GPS/GLONASS/BDS were compared to assess the performance of the multi-constellation kinematic PPP in detail.The mean square errors of these differences are shown in Table 3. From Table 3, we can see that the mean square errors of GPS were 0.045 m and 0.085 m in the east and north directions, respectively, which were much larger than those of GPS/BDS.With the GPS/BDS combination, accuracy improved by 71.11% and 41.18% over GPS only in the north and east directions, respectively.However, the positioning accuracy of GPS/GLONASS was slightly poorer than that of the GPS/BDS.The reason would be that there were more visible GPS/BDS satellites than which of GPS/GLAONSS, and the PDOP of GPS/BDS was better than that of GPS/GLAONSS) (Figure 2).The accuracies of GPS/GLONASS/BDS PPP were better than those of GPS/BDS by 30.77% and 56.00%.The standard deviations (STDs) of these differences are shown in Table 4, and the reference coordinates are the positioning results of the RTK.Smaller STDs indicate better results.As shown in Table 4, the STD of the GPS-only system was relatively large because the number of visible satellites was low (Figure 2).In contrast, the STD of the GPS/GLONASS/BDS PPP was clearly the best.

Convergence Time
To evaluate the convergence performance of multi-constellation kinematic PPP, we started from the data sets in five-minute intervals (i.e., 12 times an hour) until convergence was achieved.In this study, the position filter was considered to be converged when the positioning errors reached ±0.1 m and remained within that range.Table 5 presents the mean convergence time of RTK.Figures 5 and 6 show the kinematic PPP positioning errors of the four different processing cases.In these figures, we only present the convergence performance of PPP based on six replicates because of the space limitations.The convergence time was defined as the period from the first epoch to the converged epoch (indicated as a red line in these figures).Based on these figures and the table, we found that the GPS-only PPP required substantially more time to converge than the other processing cases.In contrast, the convergence performance of the GPS/GLONASS/BDS system was the best for all three coordinate components.Relative to that of GPS-only PPP, the convergence times of GPS/GLONASS and GPS/BDS PPP improved, particularly in the east direction.This result was very similar to that of Cai et al. [22].

Conclusions
In this study, the performance of multi-constellation kinematic PPP was assessed in terms of the positioning accuracy and convergence time.A kinematic experimental dataset was processed using a self-developed software based on the multi-constellation kinematic PPP model, Kalman filter and stochastic models.Based on the discussions above, we drew the following conclusions: The difference between the PPP kinematic and RTK solutions (GPS/BDS and GPS/GLONASS/BDS).

Conclusions
In this study, the performance of multi-constellation kinematic PPP was assessed in terms of the positioning accuracy and convergence time.A kinematic experimental dataset was processed using a self-developed software based on the multi-constellation kinematic PPP model, Kalman filter and stochastic models.Based on the discussions above, we drew the following conclusions: 1.
The availability of multi-constellation kinematic PPP was significantly better than that of GPS-only PPP because more satellites were observable and because the PDOP was better.

2.
The accuracy of multi-constellation kinematic PPP was greater than that of GPS-only, GPS/BDS and GPS/GLONASS PPP.The positioning accuracy of GPS/GLONASS PPP was slightly lower than that of GPS/BDS PPP.

3.
The convergence performance of GPS/GLONASS/BDS kinematic PPP was the best for all three coordinate components, particularly in the East direction.4.
Multi-constellation kinematic PPP can provide a positioning service with centimetre-level accuracy for dynamic users.

Figure 3 .
Figure 3. Coordinate increments in the following directions: north, east and up.

Figure 3 .
Figure 3. Coordinate increments in the following directions: north, east and up.

Figure 5 .
Figure 5.The difference between the PPP kinematic and RTK solutions (GPS and GPS/GLONASS).

Figure 6 .
Figure 6.The difference between the PPP kinematic and RTK solutions (GPS/BDS and GPS/GLONASS/BDS).
Figure 6.The difference between the PPP kinematic and RTK solutions (GPS/BDS and GPS/GLONASS/BDS).

Table 2 .
Models used in multi-constellation kinematic precise point positioning (PPP) data preprocessing.PCO: phase centre offsets; PCV: phase centre variations.

Table 3 .
The mean errors of the four processing cases (unit: m 2 ).STD: standard deviation.

Table 4 .
STDs of the positioning errors for the four processing cases (unit: m).

Table 5 .
Mean convergence time of RTK.