A Comparative Study for Performance Analysis of Kinematic Multi-Constellation GNSS PPP in Dynamic Environment

: This case study aims to investigate the e ﬀ ect of di ﬀ erent Multi-GNSS EXperiment (MGEX) precise products provided by International GNSS Service (IGS) Analysis Centers (ACs) on post-processing kinematic Precise Point Positioning (PPP) accuracy performance with di ﬀ erent satellite system combinations in a dynamic environment. Within this frame, a test was carried out in a lake and kinematic data were collected over 6 h at 1 Hz rate from the available Global Navigation Satellite System (GNSS) constellations with the geodetic-grade receiver ﬁxed on a marine vehicle for bathymetric mapping. PPP-derived coordinates were determined by a commercial GNSS post-processing software with di ﬀ erent processing approaches as GPS (Global Positioning System)-only, GPS + GLObal’naya NAvigatsionnaya Sputnikovaya Sistema (GLONASS), GPS + GLONASS + European Global Navigation Satellite System (Galileo), GPS + GLONASS + Chinese Global Navigation Satellite System (BeiDou), and GPS + GLONASS + Galileo + BeiDou. The PPP coordinates were then compared to the reference coordinates obtained from the post-processed carrier phase-based di ﬀ erential kinematic solutions. In general, the results showed that the kinematic multi-constellation GNSS PPP technique could provide positioning accuracy from cm to decimeter level as depending on the collected data constellations and used precise products in the processing. Among all solutions, the GPS + GLONASS + Galileo + BeiDou combination with German Research Centre for Geosciences (GFZ)’s precise products presented the best multi-GNSS PPP performance, rather than the other combinations and quad-constellation alternatives using di ﬀ erent precise products. In this study, the test procedure and the obtained results are given in detail.


Introduction
In recent years, the world of satellite-based positioning has experienced dramatic changes. A few decades ago, it was necessary to make carrier-phase based measurements with at least two geodetic-types of Global Navigation Satellite System (GNSS) receivers (one as a rover and at least one for the base) for achieving high accuracy positioning. However, nowadays, new methods and techniques using a single GNSS receiver have been developed to make it possible to determine the positions accurately without using data from any reference station. Precise Point Positioning (PPP), as one of the most common and successfully used methods, allows positioning with centimeter to decimeter (dm) level accuracy in static and kinematic applications by using a single GNSS receiver after convergence. Convergence (initialization time) is inevitable due to the phase ambiguity solution required to achieve a highly accurate position. However, convergence is a problem only in real-time PPP applications or where only one-way (forward or backward) post-processed PPP is applied.
Global Navigation Satellite System (BDS), and other regional satellite systems. The availability of these precise products provides multi-constellation GNSS PPP concept and allows determination of more accurate, reliable, and robust positioning by using global or regional navigation satellite systems. Within this frame, IGS as a main data provider initiated the IGS Multi-GNSS Experiment (MGEX) Project in 2012 for the tracking, collecting, and analyzing the available GNSS signals [6]. The number of institutions and multi-GNSS stations contributing to the MGEX project has been increasing over time. As of May 2020, 319 IGS multi-GNSS stations spread almost all over the world have been involved in MGEX projects [7]. The geographical distribution of these stations is given in Figure 1. Nowadays, seven IGS-ACs (Centre National d'Etudes Spatiales/Collecte Localisations Satellites (CNES/CLS), Center for Orbit Determination in Europe (CODE), German Research Centre for Geosciences (GFZ), Japan Aerospace Exploration Agency (JAXA), Shanghai Astronomical Observatory (SHAO), Technical University of Munich (TUM), and Wuhan University (WHU)) that took part in the MGEX Project have been routinely producing and distributing multi-constellation products including broadcast ephemerides (BRD), precise ephemerides (SP3) and clock products (CLK), Earth orientation parameters (ERP), IGS stations coordinates (SNX), and biases (BIA) for global or regional navigation satellite systems. The routine contributions of the above-mentioned analysis centers are given in Table 1 (as of May 2020). More information about MGEX is available in [7,9,10]. Table 1. The analysis centers of MGEX's and available products [7]. Nowadays, seven IGS-ACs (Centre National d'Etudes Spatiales/Collecte Localisations Satellites (CNES/CLS), Center for Orbit Determination in Europe (CODE), German Research Centre for Geosciences (GFZ), Japan Aerospace Exploration Agency (JAXA), Shanghai Astronomical Observatory (SHAO), Technical University of Munich (TUM), and Wuhan University (WHU)) that took part in the MGEX Project have been routinely producing and distributing multi-constellation products including broadcast ephemerides (BRD), precise ephemerides (SP3) and clock products (CLK), Earth orientation parameters (ERP), IGS stations coordinates (SNX), and biases (BIA) for global or regional navigation satellite systems. The routine contributions of the above-mentioned analysis centers are given in Table 1 (as of May 2020). More information about MGEX is available in [7,9,10].

Institution/ID Constellations Products and Intervals
When compared to single system PPP, multi-constellation GNSS PPP is expected to increase the positioning performance significantly, since it benefits from better optimized satellite (spatial) geometry and hence improves positioning availability and continuity. Thus, the shortening of the convergence time would be possible since the number of observed satellites increases in the multi-constellation GNSS PPP. This advantage is even more important in studies in which the observations are carried out in tough environments, have limited number of tracked satellites, blocked or sheltered GNSS signals, or poor satellite geometry for positioning [5,11]. Hence, it can be summarized the benefits of multi-GNSS as in the following [12,13]: • Increased number of visible satellites with stronger satellite geometry. • Improved reliability and redundancy. • Improved position and time accuracies.
• Improved positioning availability and continuity. • Improved signal-to-noise ratio. • Reduced signal acquisition time.

•
Reduced systematic biases in GNSS observations for positioning application.

•
Reduced problems, which may be caused by bulky structures such as buildings and trees. Table 1. The analysis centers of MGEX's and available products [7]. The main disadvantages of multi-constellation GNSS PPP positioning are that all the satellite systems have their own time scales and reference frames; moreover, the channels are biased and there are inter-systems biases [14]. These disadvantages can be overcome by using MGEX products.

Institution/ID
The conventional GNSS code pseudo-range and carrier-phase observation equations that are also the basic mathematical model of the GNSS-PPP algorithm are given by Equations (1) and (2). Carrier-phase and code pseudo-range observations on ith frequency for satellite j can be expressed as in [5]: where j: any GNSS satellite (G, R, E, or C), P i and Φ i : GNSS pseudorange and phase measurements on ith frequency (m), ρ: true geometric range (m), c: speed of light in vacuum (m/s), dt: receiver clock error (s), dT: satellite clock offset (s), d orb : satellite orbit error (m), d trop : tropospheric delay (m), d ion/L i : ionospheric delay on ith frequency (m), B i : phase ambiguity term on ith frequency (m), which includes the receiver and satellite phase biases and phase hardware delay biases, b P i : hardware delay bias in the code observations on ith frequency (m), and ε P i and ε Φ i : code and phase observations noises including multipath (m), respectively. On the other hand, the PPP technique uses ionosphere-free combined observables to eliminate the first-order ionospheric delay errors. In this case, after applying the MGEX precise satellite orbit and clock corrections, the ionosphere-free combined code and carrier-phase observables for a satellite j are expressed as [5]: where P IF and Φ IF : ionosphere-free combination of pseudo-range and phase measurements (m); f 1 and f 2 : two carrier-phase frequencies (Hz). The code hardware delay biases (b P i ) seen in the original observation equations do not appear in the ionosphere-free observables (as b P IF ) because they are absorbed into the receiver clock error for the GPS and into the system time difference for the other satellite constellations [5]. For an ambiguity-float PPP solution using ionosphere-free combined observations given in Equations (3) and (4), the ionosphere-free ambiguity parameter is estimated as a real-value constant.
In order to get the cm-level accurate positioning in multi-constellation GNSS PPP, the satellite and receiver antenna offsets and directional variation corrections; solid earth, ocean tide loading displacement effects; phase wind-up effects; relativistic effects; polar motion; and the others, should also be considered.
The concept of post-processing multi-constellation GNSS PPP has been handled by many researchers in static and/or kinematic applications. Within these studies, contributions of various GNSS satellite constellation combinations to the performance of the post-processing PPP utilizing with different precise products produced by many ACs were investigated with respect to positioning accuracy and convergence time [5,11,12,[15][16][17][18][19][20][21][22][23][24]. Their results indicate that multi-GNSS can significantly increase the number of observed GNSS satellites while optimizing the spatial geometry using multi-frequency data. They also concluded that, in general, additional satellite constellations improved the positioning accuracy of PPP compared to that achieved by using a single system. These experimental results also demonstrate that the multi-GNSS configuration reduced the convergence time when compared to single system PPP as an important factor limiting the PPP in some surveying applications. Furthermore, their findings imply that the precise products used during the processing of the data significantly affect the results especially in terms of accuracy. It is worth mentioning that, the datasets within these studies were generally retrieved from the permanent stations of the International GNSS Service (IGS), IGS Multi-GNSS Experiment (MGEX), EUREF Permanent Network (EPN), International GNSS Monitoring & Assessment Service (iGMAS), National Continuously Operating Reference Stations (CORS), and so on. The results given in these studies were based on the processing of the static data as kinematic PPP solutions. As a reminder, generally reference stations for IGS and similar networks are often established in ideal locations where the factors that will negatively affect measurements (low multipath and signal loss) will be minimal and that they can make observations from the all-possible satellites.
On the other hand, multi-constellation GNSS PPP in kinematic applications in a realistic scenario were investigated in a few studies. Anquela et al. [25] investigated combinations of GPS and GLONASS dual-frequency measurements for both static and kinematic (as car and walk trajectory) applications. Their results concluded that the addition of GLONASS to GPS improved satellite availability and geometry by more than 20%. Furthermore, GPS and GLONASS provided better accuracy compared to GPS-only solutions by processing magicGNSS online GNSS processing software. Choy et al. [26] made kinematic experiments with an aircraft lasting about 4 h. According to this study, adding the GLONASS constellation to the GPS improved the kinematic PPP accuracy and precision. Their results also reveal that, combined GPS and GLONASS observations shortened the convergence time compared to GPS-only. In another study, additional GLONASS constellation data improved the kinematic PPP performance when compared to GPS-only but there was no improvement with additional Galileo constellations [27]. The performance of quad-constellation PPP was investigated in [5] and they reveal that the performance of multi-GNSS PPP positioning showed better performance compared to single and dual-constellations. Yu and Gao [28] carried out kinematic measurements collected from three receivers mounted on a vehicle by observing GPS, BDS, and GLONASS constellations. After processing of the collected data using self-developed software in PPP mode, they found that triple-constellation kinematic PPP obtained better results than those of both GPS-only and dual-constellation PPP in terms of positioning accuracy and convergence time. Guo et al. [29] examined the feasibility of the multi-GNSS PPP in agriculture studies. According to processing of their measurements (static and kinematic), they concluded that multi-GNSS PPP showed better precision and repeatability than GPS-only PPP. Erol and Alkan [30] tested the performance of kinematic multi-GNSS PPP technique with different groups of constellations in a marine environment and showed that quad-constellation solution increases the position accuracy compared to other combinations. In [31] the performance of the multi-GNSS kinematic PPP in UAV photogrammetry was analyzed. The study showed that the results obtained from PPP were comparable to those obtained from the classical relative RTK method.
In this case study, unlike the above-mentioned studies, a realistic kinematic experiment was carried out on a dam lake situated in a deep valley and surrounded by high hills covered with dense trees. The performance of the kinematic multi-constellation GNSS PPP was analyzed in terms of the accuracy using different GNSS constellations (single, dual, triple, and quad-constellation) with different precise products. The multipath and signal loss due to obstruction of satellite signals were expected through the kinematic experiment that were carried out in such a challenging environment. The results of this study will guide in the selection of the most suitable satellite constellation and precise products for the more accurate PPP solution especially in inland water and marine surveys with challenging environmental conditions and severe terrain obstructions. However, as of May 2020, there are roughly 100 operational satellites in orbit from global satellite systems, and soon about 120 satellites and more will be available when all global and regional satellite systems reach full constellations, which means that any user typically can observe on average more than 20 satellites at any time. From this point of view, it is now usually possible to observe more satellites than users need and the important factor is which satellite constellations and which precise products should be used for an accurate and reliable solution. ) for Network RTK. The R10 provides ±2 cm RMS horizontal and ±5 cm RMS vertical accuracy with the Trimble CenterPoint RTX correction service as RT-PPP using satellite delivery or internet connection anywhere in the world. More information about the Trimble R10 receiver can be found in [32]. The kinematic test data were collected for approximately 6 h (about 21,350 measurement epochs) from all visible GPS, GLONASS, Galileo, and BDS satellites in view with 1-s positioning rate. The satellite elevation mask angle was set to 10 • . The trajectory of the surveying vessel during the kinematic test is shown in Figure 2, right.

The Test Area and Data Collection
vertical accuracy with the Trimble CenterPoint RTX correction service as RT-PPP using satellite delivery or internet connection anywhere in the world. More information about the Trimble R10 receiver can be found in [32]. The kinematic test data were collected for approximately 6 h (about 21,350 measurement epochs) from all visible GPS, GLONASS, Galileo, and BDS satellites in view with 1-s positioning rate. The satellite elevation mask angle was set to 10°. The trajectory of the surveying vessel during the kinematic test is shown in Figure 2, right.  The number of the tracked satellites and the average number of tracked satellites and calculated PDOP (Position Dilution of Precision) values during the kinematic test were given in Figure 3 and Table 2, respectively. According to Figure 3 and Table 2 The number of the tracked satellites and the average number of tracked satellites and calculated PDOP (Position Dilution of Precision) values during the kinematic test were given in Figure 3 and Table 2, respectively. According to Figure 3 and Table 2   In order to give an idea about the collected data quality, the Multipath and SNR (Signal-to-Noise Ratio) graphics for the carrier frequencies in the constellations were created and given as sky-plot in the Figure 4. According to plots in Figure 4, it is clearly seen that all SNR values are over the limit that is 35 dBHz, except GPS L2, and the multipath values are less than 0.60 m for both frequencies in all constellations. According to the collected data, Figure 4 shows that GPS L1 frequency signal is more powerful than L2 frequency, and L2 multipath values are bigger than L1 values. Both Galileo signals have the best SNR values and small multipath values among all constellations.  In order to give an idea about the collected data quality, the Multipath and SNR (Signal-to-Noise Ratio) graphics for the carrier frequencies in the constellations were created and given as sky-plot in the Figure 4. According to plots in Figure 4, it is clearly seen that all SNR values are over the limit that is 35 dBHz, except GPS L2, and the multipath values are less than 0.60 m for both frequencies in all constellations. According to the collected data, Figure 4 shows that GPS L1 frequency signal is more powerful than L2 frequency, and L2 multipath values are bigger than L1 values. Both Galileo signals have the best SNR values and small multipath values among all constellations. Ratio) graphics for the carrier frequencies in the constellations were created and given as sky-plot in the Figure 4. According to plots in Figure 4, it is clearly seen that all SNR values are over the limit that is 35 dBHz, except GPS L2, and the multipath values are less than 0.60 m for both frequencies in all constellations. According to the collected data, Figure 4 shows that GPS L1 frequency signal is more powerful than L2 frequency, and L2 multipath values are bigger than L1 values. Both Galileo signals have the best SNR values and small multipath values among all constellations.

Data Processes and Numerical Results
One of the goals of this case study was to investigate the effect of additional GNSS constellations on the kinematic PPP. In order to assess the accuracy performance of the PPP solutions, the collected data were processed under the following five different solution scenarios based on combining the observations from different satellite constellations with CODE's COM ID precise products: In order to obtain PPP-derived coordinates, the collected kinematic data were processed using GrafNav GNSS Post-Processing Software using the parameters given in Table 3. This software supports the GPS, GLONASS, Galileo, BDS, and QZSS global and regional satellite systems. During the kinematic PPP process, "Multi-Pass" processing strategy was chosen. In this approach, the PPPderived coordinates were calculated three times sequentially as forward, reverse, and forward again using the Extended Kalman Filter (EKF) algorithm [33]. After finalizing the processes along each direction, the converged Kalman Filter error states have been applied to the next processing direction [33]. In the EKF algorithm, 3D position, receiver clock error, tropospheric delay, and carrier-phase float ambiguities are estimated as EKF states [34]. The deterministic and stochastic models of the PPP processing and the details of the EKF algorithm used to estimate the unknowns are given in Guo and Zhang [35], Elsobeiey and El-Rabbany [36], and Zhao [37]. In the PPP module of GrafNav software, the tropospheric delay is corrected using the UNB model developed by the University of New

Data Processes and Numerical Results
One of the goals of this case study was to investigate the effect of additional GNSS constellations on the kinematic PPP. In order to assess the accuracy performance of the PPP solutions, the collected data were processed under the following five different solution scenarios based on combining the observations from different satellite constellations with CODE's COM ID precise products: In order to obtain PPP-derived coordinates, the collected kinematic data were processed using GrafNav GNSS Post-Processing Software using the parameters given in Table 3. This software supports the GPS, GLONASS, Galileo, BDS, and QZSS global and regional satellite systems. During the kinematic PPP process, "Multi-Pass" processing strategy was chosen. In this approach, the PPP-derived coordinates were calculated three times sequentially as forward, reverse, and forward again using the Extended Kalman Filter (EKF) algorithm [33]. After finalizing the processes along each direction, the converged Kalman Filter error states have been applied to the next processing direction [33].
In the EKF algorithm, 3D position, receiver clock error, tropospheric delay, and carrier-phase float ambiguities are estimated as EKF states [34]. The deterministic and stochastic models of the PPP processing and the details of the EKF algorithm used to estimate the unknowns are given in Guo and Zhang [35], Elsobeiey and El-Rabbany [36], and Zhao [37]. In the PPP module of GrafNav software, the tropospheric delay is corrected using the UNB model developed by the University of New Brunswick [38]. However, the wet part of the tropospheric delay is highly variable and cannot be modeled with sufficient accuracy. Thus, when estimating the receiver position and other unknowns, the tropospheric delay is also predicted [34]. As it was explained before, the PPP technique uses precise products in addition to collected GNSS data to estimate the accurate coordinates. Precise satellite orbit and clock corrections of CODE with COM ID retrieved from the related ftp sites were utilized to estimate the PPP-derived coordinates. Here it is noted that the CODE is one of the IGS Analysis Centers, and broadcasts the precise products for GPS, GLONASS, Galileo, BDS, and QZSS satellite constellations under the IGS MGEX Project. Table 3. Used methods and models for Precise Point Positioning (PPP) processing.

Items
Models/Methods In addition to the assessment of the accuracy performance of multi-constellation GNSS PPP, the effect of different MGEX precise products, produced by different IGS ACs on the PPP results, was also investigated. In this respect, the collected kinematic GNSS data were also processed by applying the following solution scenarios: S.6-GPS-only observations (with IGS's precise products), S.7-GPS-only observations (with CODE's precise products-COD ID), S.8-GPS+GLONASS+Galileo+BDS observations (with GFZ's precise products), and S.9-GPS+GLONASS+Galileo+BDS observations (with Wuhan University's precise products).

PPP model
The PPP-derived coordinates were also calculated by applying the aforementioned PPP solution approach with GrafNav GNSS Post-Processing Software for Scenarios 6-9. The precise GNSS satellite orbit and clock products were provided from IGS, CODE, GeoForschungsZentrum Potsdam (GFZ), and Wuhan University's Analysis Centers.
After the kinematic PPP processing given in Scenarios 1-5, the used satellites in PPP and calculated PDOP values were given in Figure 5, left and right, respectively. Figure 5 shows top-down satellite numbers (left) and PDOPs (right) for G, GR, GRC, GRE, and GRCE combinations, respectively. According to Figures 3 and 5, the number of tracked satellites and the number of satellites used in kinematic PPP processing is similar. However, the PDOP values have slightly changed. As seen from the right side of Figure 5, the average PDOP values are 2.5 for GPS, 1.8 for GPS+GLONASS, 1.4 for GPS+GLONASS+BDS, 1.4 for GPS+GLONASS+Galileo, and 1.2 for GPS+GLONASS+BDS+Galileo combinations, respectively. In the PPP processing, the change of PDOP values means that some code and phase observations or satellites are excluded from the adjustment process as an outlier because of GNSS un-modelled error sources. If the PDOP value is increasing, then the satellite number or observation number is decreasing, which means that 3D position accuracy is reduced. According to Figure 5, it is seen that every new satellite system added to the GPS system decreases the PDOP value. Under normal conditions, if PDOP values are decreasing, 3D position accuracy is expected to increase. However, different results may occur in PPP processing due to un-modeled GNSS error sources like cycle-slip and multipath. In addition, the accuracy of precise products, and the positions of satellites in the sky (distributions), used in the PPP processing, affect the accuracy of the results. The PPP-derived coordinates were also calculated by applying the aforementioned PPP solution approach with GrafNav GNSS Post-Processing Software for Scenarios 6-9. The precise GNSS satellite orbit and clock products were provided from IGS, CODE, GeoForschungsZentrum Potsdam (GFZ), and Wuhan University's Analysis Centers.
After the kinematic PPP processing given in Scenarios 1-5, the used satellites in PPP and calculated PDOP values were given in Figure 5, left and right, respectively. Figure 5 shows top-down satellite numbers (left) and PDOPs (right) for G, GR, GRC, GRE, and GRCE combinations, respectively. According to Figures 3 and 5, the number of tracked satellites and the number of satellites used in kinematic PPP processing is similar. However, the PDOP values have slightly changed. As seen from the right side of Figure 5, the average PDOP values are 2.5 for GPS, 1.8 for GPS+GLONASS, 1.4 for GPS+GLONASS+BDS, 1.4 for GPS+GLONASS+Galileo, and 1.2 for GPS+GLONASS+BDS+Galileo combinations, respectively. In the PPP processing, the change of PDOP values means that some code and phase observations or satellites are excluded from the adjustment process as an outlier because of GNSS un-modelled error sources. If the PDOP value is increasing, then the satellite number or observation number is decreasing, which means that 3D position accuracy is reduced. According to Figure 5, it is seen that every new satellite system added to the GPS system decreases the PDOP value. Under normal conditions, if PDOP values are decreasing, 3D position accuracy is expected to increase. However, different results may occur in PPP processing due to un-modeled GNSS error sources like cycle-slip and multipath. In addition, the accuracy of precise products, and the positions of satellites in the sky (distributions), used in the PPP processing, affect the accuracy of the results.

Results and Discussion
With the aim of providing an assessment on the accuracy of different PPP solutions, the reference trajectory (or known coordinates) should be obtained. In order to achieve this purpose, a nearby single reference station with precisely known ITRF coordinates was used and the coordinates of each measurement epoch were obtained with the PPK (Post-Processed Kinematic) differential method using GrafNav Software by applying "Both" processing strategy. In this approach, forward and

Results and Discussion
With the aim of providing an assessment on the accuracy of different PPP solutions, the reference trajectory (or known coordinates) should be obtained. In order to achieve this purpose, a nearby single reference station with precisely known ITRF coordinates was used and the coordinates of each measurement epoch were obtained with the PPK (Post-Processed Kinematic) differential method using GrafNav Software by applying "Both" processing strategy. In this approach, forward and reverse solutions were processed independently, and then automatically combined. This combination maximizes the solution accuracy and assists in quality control [33]. At the end of these stages, the coordinates of each measurement epoch were calculated within the precision of 15 mm and 25 mm for 2D position and height coordinates on average, respectively. The differential solution was used as reference coordinates for later comparison of the PPP solutions.
The PPP-derived coordinates obtained from all mentioned processing scenarios were compared to the known coordinates epoch      In order to visualize the distribution of the calculated differences given in Figures 6 and 7, histograms were constructed for each solution (Figures 8 and 9). In order to visualize the distribution of the calculated differences given in Figures 6 and 7, histograms were constructed for each solution (Figures 8 and 9).   Through the histograms, mean values and outliers can easily be identified in the datasets. When Figures 8 and 9 were interpreted, it was noticed that there were many outliers for the differences. If these outliers and mean values were excluded from the dataset, the differences given in Figure 8 decreased to the level of −/+10 cm for both easting and northing components, respectively. When the Figure 9. Histograms of the calculated differences for Easting (a), Northing (b) and Height (c) components, obtained from processing Scenarios 6-9.
Through the histograms, mean values and outliers can easily be identified in the datasets. When Figures 8 and 9 were interpreted, it was noticed that there were many outliers for the differences. If these outliers and mean values were excluded from the dataset, the differences given in Figure 8 decreased to the level of ±10 cm for both easting and northing components, respectively. When the differences for height coordinates were considered, it was seen that the differences decreased to the level of ±20 cm approximately (in 95% confidence interval).
When looked at Figure 9 for the processing scenarios from 6 to 9, the distributions (except G+R+E+C solution using Wuhan University's precise products) are between intervals approximately −5 cm and +5 cm for both easting and northing components and −10 cm and +10 cm for height component.
In order to make an accuracy estimate of the kinematic PPP solutions, mean values, and Root Mean Square Error (RMSE) of differences for the 2D position and height components were calculated and given in the following Tables 4-6 for all solutions. Table 4. Statistical analysis of the differences in 2D position and height components for different satellite combinations.

Scenario 5 (G+R+E+C)
using CODE's COM ID precise products 4 ±5 2 ±8 Table 5. Statistical analysis of the differences in 2D position and height components for different precise products (for GPS-only Solutions).

Scenario 7 (G)
using CODE's COD ID precise products 2 ±3 4 ±6 Table 6. Statistical analysis of the differences in 2D position and height components for different precise products (for quad-constellation Solutions).

Scenario 9 (G+R+E+C)
using Wuhan University' precise products 7 ±8 8 ±11 As given in Figure 6, Figure 8 and Table 4, the mean values reached 6 cm in the 2D position and 9 cm in height component for the GPS-only PPP processing results. The RMSE for the GPS-only are ±6 cm and ±12 cm for 2D position and height components, respectively (Scenario 1). The combination of GPS and GLONASS (Scenario 2) also produced similar results with the mean values being 6 cm in the 2D position and 11 cm in height component. When the results were examined by means of RMSE, the values obtained were ±7 cm for 2D position and ±14 cm for height components. These results imply that the solutions produce almost the same level of accuracy for both GPS-only and GPS and GLONASS constellation solutions. Because of the challenging environment, the DOP (Dilution of Precision) values of GLONASS only constellation was always over three during the kinematic test.
When the combination of the GPS+GLONASS+BDS were considered (Scenario 3), the mean values reached 4 cm in the 2D position and 1 cm in height components. The triple-constellation GNSS PPP results have the best RMSE values among them with ±5 cm for 2D position and ±7 cm for height components (Table 4). It is clearly seen that BDS constellation, in addition to GPS and GLONASS constellations, improves the 2D positioning accuracy by 28 % and 50% in 2D position and height with respect to RMSE values, respectively. When compared to the GPS-only solution, improvements of about 17% and 42% in 2D position and height, respectively, are seen.
For the other triple-constellation kinematic PPP, GPS+GLONASS+Galileo (Scenario 4), the solution provided mean values of 5 cm with ±6 cm RMSE for 2D position while 12 cm mean with ±14 cm RMSE with for height component. When the RMSE values were investigated, the 2D positional accuracy was found to be almost the same level for this solution in comparison with the both GPS-only and GPS+GLONASS combination solutions. This result indicates that the contribution of the Galileo satellites to the GPS-only and GPS+GLONASS combinations did not provide significant improvement. However, when the results obtained from GPS+GLONASS+Galileo solution were compared to GPS+GLONASS+BDS, the RMSE value was found at almost the same level for 2D position, whereas it was worse than by a factor of about two for the height component. Due to the harsh environmental conditions, the number of tracked Galileo satellites was below four in 12% of the kinematic tests (Figure 3), and the calculated DOP values were above six in 40% of the test times ( Figure 10). improvement. However, when the results obtained from GPS+GLONASS+Galileo solution were compared to GPS+GLONASS+BDS, the RMSE value was found at almost the same level for 2D position, whereas it was worse than by a factor of about two for the height component. Due to the harsh environmental conditions, the number of tracked Galileo satellites was below four in 12% of the kinematic tests (Figure 3), and the calculated DOP values were above six in 40% of the test times ( Figure 10). Finally, the quad-system, GPS+GLONASS+Galileo+BDS (Scenario 5), revealed differences between PPP to reference solution with 4 cm mean and ±5 cm RMSE for 2D position and 2 cm mean and ±8 cm RMSE for height component. This combination has almost the same mean and RMSE values with the combination of GPS+GLONASS+BDS. Looking at the RMSE values for the GPS+GLONASS+Galileo+BDS solution given in Table 4, it was seen improvements of 17%, 29%, and 17%, in the 2D position, respectively, in comparison with GPS-only, GPS+GLONASS, and GPS+GLONASS+Galileo combinations, respectively. For height component, improvements of 33%, 43%, and 43% in comparison with GPS-only, GPS+GLONASS, and GPS+GLONASS+Galileo combinations, respectively. Finally, the quad-system, GPS+GLONASS+Galileo+BDS (Scenario 5), revealed differences between PPP to reference solution with 4 cm mean and ±5 cm RMSE for 2D position and 2 cm mean and ±8 cm RMSE for height component. This combination has almost the same mean and RMSE values with the combination of GPS+GLONASS+BDS. Looking at the RMSE values for the GPS+GLONASS+Galileo+BDS solution given in Table 4, it was seen improvements of 17%, 29%, and 17%, in the 2D position, respectively, in comparison with GPS-only, GPS+GLONASS, and GPS+GLONASS+Galileo combinations, respectively. For height component, improvements of 33%, 43%, and 43% in comparison with GPS-only, GPS+GLONASS, and GPS+GLONASS+Galileo combinations, respectively.
When the differences given in Figure 7 and Table 5 were analyzed, the mean values 2 cm for 2D position and 5 cm for height were obtained from the GPS-only solution using IGS products (Scenario 6). The RMSE were calculated ±3 cm and ±6 in 2D position and height components, respectively. When the CODE's COD ID products was used (Scenario 7), the same differences were found as in the previous scenario (Scenario 6). The results imply that the positioning accuracy of GPS-only PPP using IGS precise products is comparable to using CODE analyzing center's COD-ID precise products in terms of accuracy. However, if scenarios 1 and 7 are compared for the GPS-only PPP solution, the results of two different precise products (COM-ID and COD-ID), belonging to the same institution (CODE), are two times different from each other (Table 5). Both CODE COM-ID and CODE COD-ID products are produced by the institution "Center for Orbit Determination in Europe" that is an IGS AC. COM-ID precise products include all satellite constellations and produced in the scope of the MGEX project. However, COD-ID products only support GPS and GLONASS satellite constellations. These two products group are created by using different data sets, ground stations, and strategies. Because of that, the accuracy of these products is different and they give different PPP results for the same data and constellation.
On the other hand, according to the results obtained from the quad-constellation system using the products of the GFZ (Scenario 8), mean difference of 3 cm was obtained for the 2D position and 0 cm for the height component. Using the Wuhan University products (Scenario 9), the mean values were obtained 7 cm for the 2D position and 8 cm for height. According to Figure 9 and Table 6, the RMSE values were obtained as ±4 cm for 2D position and ±5 cm height component for Scenario 8, and ±8 cm for 2D position and ±11 cm height components for Scenario 9, respectively. This result implies that, the solutions produce different levels of accuracy depending on the used different precise products.
According to comparisons of the quad-constellation used with CODE'S COM products (Scenario 5), GFZ's precise products (Scenario 8), and Wuhan University's precise products (Scenario 9), the best results were obtained from GFZ products when the mean values between PPP solution and reference coordinates were considered. In terms of accuracy, GPS+GLONASS+Galileo+BDS with GFZ's precise products also present the best performance rather than the other quad-constellation alternatives. It is seen in Tables 4-6 that the quad-system using GFZ's precise products achieves a better PPP solution in terms of mean value and accuracy (RMSE) when compared to the other constellation combinations and precise products provided by different ACs.
The epoch-by-epoch 2D position and ellipsoidal height differences between the PPP and PPK derived coordinates for the best multi-GNSS PPP solution (Scenario 8) are given in the Figure 11 as a 2D plot on the surveying path. were obtained 7 cm for the 2D position and 8 cm for height. According to Figure 9 and Table 6, the RMSE values were obtained as ±4 cm for 2D position and ±5 cm height component for Scenario 8, and ±8 cm for 2D position and ±11 cm height components for Scenario 9, respectively. This result implies that, the solutions produce different levels of accuracy depending on the used different precise products.
According to comparisons of the quad-constellation used with CODE'S COM products (Scenario 5), GFZ's precise products (Scenario 8), and Wuhan University's precise products (Scenario 9), the best results were obtained from GFZ products when the mean values between PPP solution and reference coordinates were considered. In terms of accuracy, GPS+GLONASS+Galileo+BDS with GFZ's precise products also present the best performance rather than the other quad-constellation alternatives. It is seen in Tables 4-6 that the quad-system using GFZ's precise products achieves a better PPP solution in terms of mean value and accuracy (RMSE) when compared to the other constellation combinations and precise products provided by different ACs.
The epoch-by-epoch 2D position and ellipsoidal height differences between the PPP and PPK derived coordinates for the best multi-GNSS PPP solution (Scenario 8) are given in the Figure 11 as a 2D plot on the surveying path.  Considering the histograms of all calculated differences, the differences of CODE's COM ID PPP fell within the range of approximately −10 cm and +10 cm for both easting and northing components, respectively, and approximately −20 cm to +20 cm for height component ( Figure 8). However, the differences in PPP solutions using other IGS Analysis Centers (ACs) precise products ( Figure 9) are 2 times better than the previous solutions.

-
Additional constellations did not improve kinematic PPP solution results in all cases. For instance, although triple-constellation (GPS+GLONASS+BDS) PPP using CODE's COM ID provides a significant improvement over GPS-only and GPS+GLONASS, when the other triple-constellation (GPS+GLONASS+Galileo combination) was considered, no significant improvement has been achieved.

-
Similarly, quad-constellation (GPS+GLONASS+Galileo+BDS) PPP with GFZ's precise products provides the best results in all quad PPP solutions (Table 6). From this point of view, there is no doubt that additional satellites generally have improved the positioning performance; however, based on these results it is concluded that the PPP positioning accuracy greatly depends on the precise satellite orbit and clock products applied by the user.

-
As described above, the PPP coordinates were calculated using GrafNav software multi-pass strategy. In this approach, the coordinates are calculated as forward, then reverse, and then forward again. Therefore, in this study, the effect of additional GNSS constellations on the PPP convergence time was not investigated.

-
The findings imply that kinematic multi-constellation GNSS PPP technique allows to reach cm-to decimeter-level of accuracy in a challenging dynamic environment.

Conclusions
In this case study, the contribution of the additional satellite systems on kinematic PPP solution performance in terms of accuracy was investigated. In addition, the effects of different MGEX precise satellite orbit and clock products from different IGS Analyzing Centers on the PPP results were also examined. In order to assess these issues, a real kinematic test was carried out and data were collected for 6-h with 1 Hz from a GNSS receiver setup on a vessel. According to Figures 3 and 5, it can be concluded that the additional GNSS constellations improved the observed number of satellites in comparison with that of any single-system constellation. The study reveals that in general, additional GNSS constellations improve the solution accuracy by observing more satellites compared to GPS alone. This is more important when the surveying was carried out in the field where satellite signals are partially obstructed or have multipath potential or have signal loss as a result of obstruction etc.
One of the main conclusions of the kinematic test showed that the obtained accuracy was affected by not only the observed GNSS constellations but also highly dependent on the chosen precise orbit and satellite clock products that were introduced in the processing stage. As is known, there are different products provided by a variety of analysis centers, and these centers generally have different datasets, processing strategies, and software for producing their products, and thus, the accuracy and quality of the orbit and clock corrections may differ from each other. According to all solutions conducted here, the GPS+GLONASS+Galileo+BeiDou combination with GFZ's precise products presented the best multi-GNSS PPP performance rather than the other combinations and quad-constellation alternatives using different precise products.
The overall results obtained from nine different solutions show that the multi-constellation GNSS kinematic PPP can provide cm-to dm-level of 3D positioning accuracy especially in applications carried out challenging inland waters and marine conditions. Offshore positioning, precise marine engineering, construction, dredging, cable and pipeline laying, ocean buoy positioning, tsunami warning, precise hydrographic surveying, and coastal mapping are examples of these kinds of applications. In addition, the accuracies obtained from this study meet the minimum positional accuracy requirement for preparing nautical charts according to many international hydrographic organizations and authorities regulations, like the International Hydrographic Organization (IHO).
This level of accuracy also meets the requirements of most geodetic surveying applications; like dynamic vehicle positioning, autonomous vehicle systems positioning, UAV remote sensing, UAV photogrammetry, mobile mapping systems, aviation, precise agriculture, tectonic geodesy, deformation monitoring, atmospheric water vapor sensing, and low orbiter satellite positioning.
This study shows that PPP is a strong alternative to the conventional carrier phase based differential technique in terms of accuracy, ease of use, and operating cost. The study reveals that the performance of multi-constellation GNSS PPP would be improved in the near future after adding more satellites to the existing systems in orbit as well as launching of new satellite systems (such as Korean Positioning System-KPS) and setting up more multi-GNSS reference stations.
Funding: This research received no external funding.