A Robust Algorithm for Multi-GNSS Precise Positioning and Performance Analysis in Urban Environments

: In this paper, we propose a partial ambiguity method of global navigation satellite system (GNSS) reliable positioning based on a robust estimation model to address the problems of the low reliability and availability of GNSS positioning in urban complex environments. First, the high-precision observations selected on the basis of the signal-to-noise ratio (SNR) were used to solve ambiguities. Then, the ﬁxed ambiguities were used as constraints to solve the ambiguities of low-quality observations. The robust estimation method was used to reduce the impact of outliers for the ambiguity solutions. The robust estimation was also used to solve the position parameters to reduce the inﬂuence of the residual errors and uncorrected ambiguities for GNSS high-accuracy positioning. Static and dynamic data were used to evaluate the proposed algorithm. These experiments show that the proposed algorithm with the robust estimation can reduce the ﬁxed time of ambiguity initialization, compared with the conventional algorithm without the robust estimation. The positioning accuracy and solution rate are similar regardless of whether the robust estimation is used in the GNSS unblocked environment. In blocked environments, the solution rate improves to more than 99%, and the three-dimensional (3D) position accuracy improves by more than 70% when the robust estimation is used. When the observation number of simulated small gross error accounts for 40.91% of total observations, the centimeter-level positioning accuracy can still be obtained via several robust estimation models. In the urban blocked environment, the IGG (Institute of Geodesy and Geophysics) III scheme has a better performance than other robust schemes discussed in this paper with regard to the positioning performance and computational efﬁciency.


Introduction
The global network completion of the Beidou Navigation Satellite System (BDS), along with the Global Positioning System (GPS) and Global Navigation Satellite System (GLONASS), ensures that there are enough satellites currently operating [1]. The geometric structure of a single system can be improved by a multi-satellite system, which ensures positioning accuracy and reliability [2]. Multisystem global navigation satellite system (GNSS) positioning is a research trend, and it has been widely applied in mobile mapping, intelligent driving, and smart agriculture [3,4].
A high-accuracy, continuous, and robust position reference is required for mapping of mobile mapping systems and intelligent driving. However, GNSS observations are susceptible to the influence of blocked environments, and gross errors inevitably occur in observations, resulting in biased estimations [5,6]. Currently, these processing methods of gross errors are outlier detection and robust estimation [7]. The original outlier detection theory was proposed in Baarda's study, which makes a decision between the null hypothesis and alternative hypothesis [8]. The observation of outlier that was detected was excluded,

Observation Model
The double-difference method can eliminate or reduce most errors of GNSS observations, and it is widely used for high-accuracy GNSS positioning. The double-difference equations of the GNSS based on observations of the pseudo-range and carrier phase are as follows [30,31]: ∇∆P sr i,mn = ∇∆ρ sr mn + ∇∆T sr mn + µ i ∇∆I sr mn + ∇∆d sr mn + ∇∆ε sr i,mn,ρ , λ i ∇∆ϕ sr i,mn = ∇∆ρ sr i,mn + ∇∆T sr mn − µ i ∆I sr mn + λ j ∇∆N sr i,mn + ∇∆d sr mn + ∇∆ε sr i,mn,ϕ , (2) where s and r represent satellite numbers, i, m, and n represent the satellite frequency number, reference station, and mobile station, respectively, ∇∆ is the double-difference operator, P represents the pseudo-range observation, ϕ represents the carrier-phase observation, ρ is the geometric distance from the satellite to the receiver, λ i is the wavelength of the ith frequency, N is the integer ambiguity, T represents the tropospheric delay, I represents the ionospheric delay of the first frequency, µ i = f 2 1 / f 2 i is the conversion coefficient of ionospheric delay from the first frequency to ith frequency, f 1 is the frequency of the first waveband, d represents other errors that are independent of the frequency, and ε ρ and ε ϕ represent noises of the pseudo-range and carrier phase, respectively.
After linearization of Equations (1) and (2), the error equation can be obtained as follows: where V represents the residual, B and C are design matrices of the parameter estimations, X represents the estimated position parameters, and L represents the observation value.
The least squares method is used to solve Equation (3), and the floated ambiguities and the corresponding variance-covariance are obtained. In this paper, we used the LAMBDA algorithm to fix ambiguities. The ratio value and unit weight of the standard error were used to confirm ambiguity. To further improve the accuracy of the ambiguity solution, we propose a partial ambiguity method based on wide-lane (WL) guidance. First, the WL ambiguity of high-quality observations selected on the basis of the signal-to-noise ratio (SNR) is fixed, and then the fixed ambiguity is used as a constraint to solve the residual ambiguities [32]: where f ix is the fixed solution, and f loat represents the floated solution.
The least squares method is used to solve Equation (4), and the strength of the equation is enhanced because of the constraint of the fixed ambiguities. The robust estimation is used in solving Equation (4) to reduce the influence of outliers for the ambiguity solution.
After the WL ambiguity is fixed, the basic ambiguity can be solved according to the relationship between the WL observation equation and the basic observation equation. The fixed WL ambiguity is brought into the WL observation equation of the carrier phase, and then the WL observation equation subtracts the basic observation to calculate basic ambiguity.
where N b, f loat is the float ambiguity of the basic observation equation, N wl, f ix is the fixed WL ambiguity, ϕ wl is the WL observation of the carrier phase, and λ wl is the WL wavelength. The float ambiguity calculated using Equation (5) can be rounded directly for a short baseline. After the basic ambiguity is fixed, it is substituted into Equation (3), and the observation equation based on the carrier phase is used to solve the position parameters: where P ϕ is the weight matrix. The SNR can reflect the quality of observations. When the quality of observations is lower, the SNR is smaller. Therefore, the random model based on the SNR was used in this paper: where S is the SNR, B represents the phase tracking bandwidth, and C i takes an empirical value in practical applications [33].
To resist outliers of the GNSS observations, the robust estimation is used to improve the accuracy of the parameter estimation: where P ϕ,i = γP i is the equivalent weight, and γ is the weight factor of the robust estimation. In urban blocked environments, different weight factors are selected, and the robust performance of GNSS positioning is different. Next, several equivalent weight functions are discussed.

Andrews Function
The Andrews function contains the weight reduction segment and elimination segment, and its weight factor can be expressed as follows [34]: where v i = v i /σ i is the standardized residual, σ i is the standard deviation, and c is a constant (taken as 1.339 in this paper).

Danish Function
The Danish function includes a normal segment and weight reduction segment, and its weight factor can be expressed as follows [7]: where exp(•) represents the exponential function, and c is a constant (taken as 1.5 in this paper).

Hampel Function
The Hampel function includes the normal segment, weight reduction segment, and eliminated segment, and its weight factor can be expressed as follows [16]: where a, b, and c are the modulation coefficients, which are 1.7, 3.4, and 8.5, respectively, in this paper.

Huber Function
The Huber function consists of a normal segment and weight reduction segment, and its weight factor can be expressed as follows [35,36]: where c is a constant (taken as 2.0 in this paper).

IGGIII Function
The IGGIII function contains the normal segment, weight reduction segment, and eliminated segment, and its weight factor can be expressed as follows [13]: where k 0 represents the quantile parameter, and k 1 represents the elimination points. In this paper, k 0 and k 1 were taken as 1.2 and 2.8, respectively.

Tukey Function
The Tukey function contains the weight reduction segment and elimination segment, and its weight factor can be expressed as follows [37]: where c is the regression factor (taken as 4.685 in this paper). Figure 1 shows the weight factor of different robust estimation functions described in Section 2.2.1 to Section 2.2.6, with v i ranging from −10 m to 10 m. In Figure 1, the curves of weight factor of Andrews function and Tukey function are similar. The parameters estimated using Equations (9) and (10) are either decreased in weight as is or the outliers are discarded, and the weight factor of observation is one only when v i is equal to zero, which almost reduces the weight of observations used for parameters estimation. The weight value calculated by the Danish function is either left as is or decreased, and the weight factor gradually approaches zero as |v i | increases. The weight factor calculated by the Hampel function includes four segments. The weight factor is one when the |v i | is less than or equal to 1.7 m, the weight factor is zero when |v i | is greater than 8.5 m, and the weight factor is decreased when |v i | is greater than 1.7 m but less than or equal to 8.5. The weight value calculated by the Huber function is either left as is or decreased, and the weight factor is still large even when |v i | is 10 m. That is, the observation with large outliers still has a certain weight value that decreases the accuracy of parameter estimation. The weight factor of the IGGIII function includes the normal segment, weight reduction segment, and eliminated segment. The weight factor is one when |v i | is less than 1.2 m, the weight factor is gradually decreased as the |v i | decreases, and the weight factor is zero when |v i | is greater than or equal to 2.8 m. Compared with other robust estimation function, the weight factor of IGGIII is smaller when |v i | is greater than 1.3.  The equivalent weight function discussed in Section 2.2 is substituted in (8) for iterative calculation to solve the position parameters [38]. The equivalent weight function discussed in Section 2.2 is substituted into Equation (8) for iterative calculation to solve the position parameters [38].
where k represents the kth iteration, and r represents the number of redundant observations. In this paper, when the difference between the position parameters calculated in the two consecutive iterations meets the convergence accuracy, the calculation is stopped.
The flowchart of the proposed algorithm is shown in Figure 2. The process can be summarized as follows: the data of base station and mobile station are used to estimate parameters. The high-quality observations are selected on the basis of the SNR to build the subset of the WL ambiguity, and the ratio value and unit weight of standard error are used to confirm the ambiguity. The WL equations of low-quality observations are constructed by the constraint of the fixed WL ambiguity, and the robust estimation is used to solve the residual WL ambiguities. The LAMBDA algorithm is used to search integer ambiguities, and the ratio value and unit weight of standard error are used to confirm the integer ambiguity. The basic ambiguities are solved according to the relationship between the wide-lane observation equation and the basic observation equation, and the fixed basic ambiguities are substituted into the observation equation to solve the position parameters. The robust estimation is used to reduce the effects of outliers for position parameter estimation. used to solve the residual WL ambiguities. The LAMBDA algorithm is used to search integer ambiguities, and the ratio value and unit weight of standard error are used to confirm the integer ambiguity. The basic ambiguities are solved according to the relationship between the wide-lane observation equation and the basic observation equation, and the fixed basic ambiguities are substituted into the observation equation to solve the position parameters. The robust estimation is used to reduce the effects of outliers for position parameter estimation.

Results
A NovAtel ProPak6 receiver was used to collect static and dynamic data to analyze the performance of the proposed algorithm. Seven schemes were designed, and selfdevelopment programs were used to analyze these data. Scheme 1 was a robust estimation based on the Andrews function. Scheme 2 was a robust estimation based on the Danish function. Scheme 3 was a robust estimation based on the Hampel function. Scheme 4 was a robust estimation based on the Huber function. Scheme 5 was a robust estimation based on the IGGIII function. Scheme 6 was a robust estimation based on the Tukey function. Scheme 7 was the method without robust estimation. The combined results of commercial software Inertial Explorer (IE) were used as a reference.

Case One
Static data were collected in an open environment. The reference station was set up on the roof, and the mobile station was set up on the square. The sample rate of reference station and mobile station was 20 Hz, and the collection time was about 3.0 h. These receivers collected data from GPS, GLONASS, and BDS. The data collection environment is shown in Figure 3.

Case One
Static data were collected in an open environment. The reference station was set up on the roof, and the mobile station was set up on the square. The sample rate of reference station and mobile station was 20 Hz, and the collection time was about 3.0 h. These receivers collected data from GPS, GLONASS, and BDS. The data collection environment is shown in Figure 3.   Figure 4, the number of visible satellites was stable during the data collection period, with more than five GPS (G) and GLONASS (R) satellites and more than nine BDS (C) satellites. The number of GPS/GLONASS/BDS (G/R/C) satellites was significantly higher than that of a single-satellite system, with the average number of G/R/C satellites increasing to more than 20.   Figure 4, the number of visible satellites was stable during the data collection period, with more than five GPS (G) and GLONASS (R) satellites and more than nine BDS (C) satellites. The number of GPS/GLONASS/BDS (G/R/C) satellites was significantly higher than that of a singlesatellite system, with the average number of G/R/C satellites increasing to more than 20.  The static data were solved by self-development programs, and the combined results of GNSS RTK of IE were used as a reference. Table 1 shows the solution rate of the seven schemes for the static data. The solution rate was the ratio of the solved epoch number of the proposed schemes to the actual epoch number. As shown in Table 1, the solution rate of the seven schemes was the same, reaching 100% regardless of whether the robust estimation was used, due to the good quality of observations.   The static data were solved by self-development programs, and the combined results of GNSS RTK of IE were used as a reference. Table 1 shows the solution rate of the seven schemes for the static data. The solution rate was the ratio of the solved epoch number of the proposed schemes to the actual epoch number. As shown in Table 1, the solution rate of the seven schemes was the same, reaching 100% regardless of whether the robust estimation was used, due to the good quality of observations.  Figure 5 shows the RMS of the position errors of several proposed schemes in the N, E, and U directions. As shown in Figure 4, the RMS of Scheme 7 (without the robust estimation) could achieve centimeter-level accuracy, at 0.0052 m, 0.0059 m, and 0.0108 m in the N, E, and U directions, respectively. The position accuracy of the robust schemes was basically the same, which was equivalent to Scheme 7 in the horizontal direction and slightly worse than Scheme 7 in the U direction. This is because the high-quality observations were also downweighted or were eliminated when these robust schemes were used in some epochs, which reduced the effect of high-quality observations for the position solution.  Table 2 shows the epoch numbers of the ambiguity initialization of several schemes for the static data. As shown in Table 2, Scheme 7 required 186 epochs to fix the ambiguity, while the robust schemes only required a single epoch to fix the ambiguity. The fixed time of ambiguity initialization was improved by more than 99%. Therefore, the robust estimation method could effectively shorten the fixed time of the ambiguity. Table 2. Epoch numbers of the ambiguity initialization of several schemes for the static data.

Scheme 1 Scheme 2 Scheme 3 Scheme 4 Scheme 5 Scheme 6 Scheme 7
Numbers of epochs To analyze the performance of different robust schemes to resist outliers, small and large gross errors were simulated in observations. Table 3 shows the epoch and size of the simulated gross errors. The outlier numbers of 1, 3, 6, 9, and 12 were simulated in observations. Specifically, one outlier was simulated in observations of GPS, and one outlier, two outliers, three outliers, and four outliers were simulated in observations of GPS, GLONASS and BDS, respectively. Figure 6 shows the position errors of different schemes of simulated gross error with different numbers and sizes. Figure 6a presents the position errors of Epoch 179,600 where gross error numbers of 0, 1, 3, 6, 9, and 12 were simulated. In this epoch, there were 22 satellites. Figure 6a shows that the GNSS positioning results were not affected when the number of simulated gross error was less than 4, and the positioning error of  Table 2 shows the epoch numbers of the ambiguity initialization of several schemes for the static data. As shown in Table 2, Scheme 7 required 186 epochs to fix the ambiguity, while the robust schemes only required a single epoch to fix the ambiguity. The fixed time of ambiguity initialization was improved by more than 99%. Therefore, the robust estimation method could effectively shorten the fixed time of the ambiguity. To analyze the performance of different robust schemes to resist outliers, small and large gross errors were simulated in observations. Table 3 shows the epoch and size of the simulated gross errors. The outlier numbers of 1, 3, 6, 9, and 12 were simulated in observations. Specifically, one outlier was simulated in observations of GPS, and one outlier, two outliers, three outliers, and four outliers were simulated in observations of GPS, GLONASS and BDS, respectively. Figure 6 shows the position errors of different schemes of simulated gross error with different numbers and sizes. Figure 6a presents the position errors of Epoch 179,600 where gross error numbers of 0, 1, 3, 6, 9, and 12 were simulated. In this epoch, there were 22 satellites. Figure 6a shows that the GNSS positioning results were not affected when the number of simulated gross error was less than 4, and the positioning error of Scheme 7 increased more than twofold when the observation number of simulated gross error accounted for 27.27% of the total observations. With the increase in the number of gross errors, the position error gradually increased. When the number of the simulated gross error was 9, these proposed robust schemes could resist gross errors. When the number of simulated gross error reached 12, these robust schemes could not effectively deal with outliers. Figure 6b presents the position errors of different schemes of Epoch 180,200. In this epoch, there were 22 satellites. As shown in Figure 6b, the positioning results of Scheme 7 were biased when the number of simulated gross errors was 1. Schemes 1, 2, 5, and 6 could still effectively resist gross errors when the number of simulated gross errors was 9. However, the positioning result was biased when the number of simulated gross errors was 9 for Schemes 3 and 4, and the position errors were abnormal when the number of simulated gross errors was 6 for Scheme 4. Figure 6c presents the position errors of different schemes of Epoch 181,400. In this epoch, the satellite number was 21. As shown in Figure 6c, the positioning results of the robust schemes (except for Scheme 1) had decimeter-level deviation; when the number of simulated gross errors was 9, Scheme 5 had the least position bias. Figure 6d presents the position errors of different schemes of Epoch 182,000. In this epoch, there were 19 satellites. As shown in Figure 6d, the positioning results of Scheme 3 and Scheme 4 were biased when the number of simulated gross errors was 6. When the observation number of simulated gross errors accounted for 47.37% of the total observations, none of the robust schemes could resist gross errors. Figure 6e presents the position errors of different schemes of Epoch 182,600. In this epoch, there were 21 satellites. As shown in Figure 6e, when the number of simulated gross errors was 6, the robust schemes (except Scheme 3 and Scheme 4) could effectively resist gross errors. When the number of simulated gross errors reached 9, these proposed robust schemes could not resist gross errors. Figure 6f presents the position errors of different schemes of Epoch 183,200. In this epoch, there were 23 satellites. In Figure 6f, when the number of simulated gross errors was 6, the robust schemes (except Scheme 4) could effectively resist gross errors. When the number of simulated gross errors was 9, robust schemes (except Scheme 1) could not resist gross errors.
In conclusion, the scheme without robust estimation could still obtain centimeter-level results when the observation number of the simulated 0.5-cycle gross error accounted for less than 13.64% of the total observations. Several robust schemes could effectively resist gross errors when the observation number of simulated gross errors accounted for 40.91% of the total observations. The robust performance of Scheme 4 was reduced when the observation numbers of the simulated one-cycle gross error accounted for 27.27% of the total observations. The robust performance of Schemes 3 and 4 was reduced when the observation numbers of simulated gross errors accounted for 40.91% of the total observations, but the position error was still at the centimeter level. Several robust schemes could still effectively resist gross errors when the observation number of the simulated five-cycle gross error accounted for 28.57%. For the simulated large gross errors of more than 10 cycles, the robust performance of Scheme 3 and Scheme 4 was slightly worse, the performance of other schemes was similar, and Scheme 1 was slightly better. Therefore, the robust effect of Schemes 3 and 4 was slightly worse for the simulated multiple gross errors. rors reached 9, these proposed robust schemes could not resist gross errors. Figure 6f presents the position errors of different schemes of Epoch 183,200. In this epoch, there were 23 satellites. In Figure 6f, when the number of simulated gross errors was 6, the robust schemes (except Scheme 4) could effectively resist gross errors. When the number of simulated gross errors was 9, robust schemes (except Scheme 1) could not resist gross errors. In conclusion, the scheme without robust estimation could still obtain centimeter-level results when the observation number of the simulated 0.5-cycle gross error accounted for less than 13.64% of the total observations. Several robust schemes could effectively resist gross errors when the observation number of simulated gross errors accounted for 40.91% of the total observations. The robust performance of Scheme 4 was reduced when the observation numbers of the simulated one-cycle gross error accounted for 27.27% of the total observations. The robust performance of Schemes 3 and 4 was reduced when the observation numbers of simulated gross errors accounted for 40.91% of the total observations, but the position error was still at the centimeter level. Several

Case Two
In urban environments, dynamic data were collected by the vehicle-born surveying system researched and development by Shandong University of Science and Technology and the company of Supersurs. The vehicle system is shown in Figure 7, and it was equipped with multiple sensors. The reference station was set up on the roof, and the mobile station was carried on the vehicle. The data included those from GPS, GLONASS, and BDS. The sample rate was 20 Hz, and the collection time was about 2.8 h. There were many blocked environments, including tall buildings, viaducts, tall porches, and shady roads. The experimental area and mobile trajectory are shown in Figure 8. In Figure 8, the marked A is the experimental area of signal interference, B is the experimental area of densely built, C is the experimental area of viaduct, and D is the experimental area of tall porch area.   Figure 9 shows the number of visible satellites. In Figure 9, the satellite num GPS was more than five, the satellite number of GLONASS was more than four, a satellite number of BDS was more than eight most of the time. The num GPS/GLONASS/BDS satellites was significantly higher than that of the single-s system, and the number of GPS/GLONASS/BDS satellites reached 15 most of the However, the satellite number of GNSS was less than four in some blocked en ments.   Figure 9 shows the number of visible satellites. In Figure 9, the satell GPS was more than five, the satellite number of GLONASS was more than satellite number of BDS was more than eight most of the time. Th GPS/GLONASS/BDS satellites was significantly higher than that of the s  Figure 9 shows the number of visible satellites. In Figure 9, the satellite number of GPS was more than five, the satellite number of GLONASS was more than four, and the satellite number of BDS was more than eight most of the time. The number of GPS/GLONASS/BDS satellites was significantly higher than that of the single-satellite system, and the number of GPS/GLONASS/BDS satellites reached 15 most of the time. However, the satellite number of GNSS was less than four in some blocked environments.
The self-development programs were used to solve the data, and the combined smoothed results of GNSS RTK/INS integrated navigation of IE were used as a reference. Table 4 shows the solution rate of the seven schemes for the dynamic data. As shown in Table 4, the solution rate of Scheme 7 (where robust estimation was not used) was 97.83%, while the solution rate of several robust schemes was improved to more than 99%. The solution rate of Scheme 6 reached 99.57%, and the improvement effect of the solution rate was the best. This was followed by Scheme 5 with a solution rate of 99.47%, while the solution rate of Scheme 1 also reached 99.35%. The self-development programs were used to solve the data, a smoothed results of GNSS RTK/INS integrated navigation of IE were u Table 4 shows the solution rate of the seven schemes for the dynamic Table 4, the solution rate of Scheme 7 (where robust estimation w 97.83%, while the solution rate of several robust schemes was impro 99%. The solution rate of Scheme 6 reached 99.57%, and the improve solution rate was the best. This was followed by Scheme 5 with a solut while the solution rate of Scheme 1 also reached 99.35%.    Table 5 shows the epoch numbers of ambiguity initialization of several schemes for the dynamic data. As shown in Table 5, 48 epochs were needed to fix the ambiguity for Scheme 7, while the ambiguity could be fixed by a single epoch when robust estimation was used. The fixed time of ambiguity initialization was improved by more than 97.92%. Thus, the robust estimation method could effectively shorten the fixed time of the ambiguity initialization.  Figure 10 shows the position error of these proposed schemes. In Figure 10a, the area marked A is the signal interference area, B is the densely built area, C is the viaduct area, and D is the tall porch area. The actual environments of these marked areas are shown in Figure 8. Figure 11 shows the 3D position errors of the marked areas in Figure 10a. Figures 10 and 11a show that Scheme 7 could not obtain the positioning results during a period of time in the area marked A, and the positioning results had decimeter-level errors in A. The positioning accuracy and solution rate were improved by robust schemes, and there were still decimeter-level errors for all robust Schemes in A. Figure 9 shows that the number of satellites during this period of time was low. The weight functions of Schemes 1 and 6 included a weight reduction segment and elimination segment. The low-quality satellites were directly eliminated such that the number of satellites participating in the position calculation was low and the reliability of positioning results was reduced. Figures  10 and 11b show that Scheme 7 had continuous decimeter errors in the area marked B, and the position errors were reduced to the centimeter level by robust schemes. As shown in Figures 10 and 11c, Scheme 7 had decimeter errors, and the position could not be obtained at some moments in the area marked C. The position accuracy was greatly improved when robust estimation was used, and Schemes 4 and 5 had a better improvement effect. As shown in Figures 10 and 11d, Scheme 7 failed to obtain the position results for a long time in the area marked D, and the solution rate was improved by the robust schemes. Due to the serious blocked environment, several robust schemes also had decimeter-level errors; in Schemes 1 and 6, the decimeter-level errors were large, whereas, in Schemes 3 and 5, the decimeter errors were small.  Table 6 shows the accuracy statistics of the position errors of the proposed schemes shown in Figure 10. As shown in Table 6 Table 6 shows the accuracy statistics of the position errors of the proposed schemes shown in Figure 10. As shown in Table 6, the RMS of the position errors of Scheme 7 was 0.0203 m, 0.0279 m, and 0.1685 m in the N, E, and U directions, respectively, and the maximum 3D position error reached 0.8113 m. The position accuracy was significantly improved by robust schemes. The RMS values of the 3D position error for the six robust schemes were 0.0516 m, 0.0493 m, 0.0495 m, 0.0496 m, 0.0494 m, and 0.0516 m. Compared with Scheme 7, the RMS of the six robust schemes was improved by 70.00%, 71.34%, 71.22%, 71.16%, 71.28%, and 70.26%, respectively. The maximum 3D position error was improved by these robust schemes. Scheme 5 had the best improvement effect, with the maximum 3D position error reduced to 0.2544 m, while the improvement effect of Schemes 1 and 6 was slightly worse.  Figure 12 shows the cumulative probability of the position errors of the proposed schemes. As shown in Figure 12, the probability of Scheme 7 was only 63.18% when the 3D position error was less than 0.05 m. The cumulative probability was improved by the robust schemes. The cumulative probability of Scheme 5 reached 71.83%, of Scheme 2 reached 71.68%, and of other robust schemes reached more than 70%. The cumulative probability of Scheme 7 was only 86.98% when the 3D position error was less than 0.10 m, and the probability of several robust schemes reached more than 94%, among which Schemes 2 and 4 reached 94.84%, and Scheme 5 reached 94.64%. The cumulative probability of Scheme 7 was only 90.51% when the 3D position error was less than 0.15 m, and the probability of several robust schemes reached more than 98%, among which Scheme 5 had the best effect with a probability of 99.15%, followed by Scheme 6 at 99.13%. Compared with Scheme 7, the cumulative probability of the position error was also improved by several robust schemes for other position error intervals.  Table 7 shows the average number of iterations of the proposed robust schemes. In this paper, the iteration calculation was terminated when the position difference of two consecutive iterations was less than a certain threshold. As shown in Table 7, the average number of iterations of Schemes 1 and 6 was more than 200, and the average number of iterations of Schemes 3 and 4 was also more than 30. The average number of Schemes 2  Table 7 shows the average number of iterations of the proposed robust schemes. In this paper, the iteration calculation was terminated when the position difference of two consecutive iterations was less than a certain threshold. As shown in Table 7, the average number of iterations of Schemes 1 and 6 was more than 200, and the average number of iterations of Schemes 3 and 4 was also more than 30. The average number of Schemes 2 and 5 was less, being only 11.3081 for Scheme 5. Therefore, Scheme 5 had the smallest computational burden among the proposed robust estimation schemes in this paper.

Discussion
In this study, the proposed strategy of GNSS positioning could obtain reliable and highaccuracy position parameters in urban blocked environments. Compared with methods without robust estimation, the performance of GNSS positioning was improved by robust schemes. Scheme 5 had the best performance among the proposed robust schemes with regard to the positioning performance and computational efficiency. The positioning performance of the proposed algorithm was analyzed on the basis of data collected in urban environments, and decimeter-level or centimeter-level positioning results could be obtained by the proposed algorithm when gross errors occurred in observations. Centimeter-level positioning accuracy could still be obtained by the proposed strategy when observation numbers of the simulated 15-cycle gross error accounted for 28.75% of the total observations. Therefore, a continuous and high-accuracy position reference can be provided for mapping by the proposed algorithm, which can promote the development of mobile mapping systems and intelligent driving in urban environments.

Conclusions
In a complex urban environment, the reliability of GNSS positioning has numerous requirements, and GNSS signals are susceptible to the interference from block environments, resulting in a reduction in the observation quality, which affects the accuracy and reliability of GNSS positioning. Therefore, we proposed an algorithm of GNSS reliable positioning with several robust schemes. The technology of double-difference GNSS RTK was used to discuss the proposed algorithm, and the least squares method was used to solve parameters. The partial ambiguities method was used to solve the ambiguities. Different equivalent weight functions of robust estimation were analyzed, and the designed schemes of robust estimation were used to reduce the influence of the residual errors and uncorrected ambiguities for high-accuracy positioning. Static and dynamic data were used to assess the performance of the proposed algorithm. Some conclusions can be drawn from our experiments: Firstly, reliable results of GNSS positioning can be obtained regardless of whether robust estimation is used in the open environment, and the accuracy and solution rate can be improved by the robust scheme in blocked environments.
Secondly, the fixed time of the ambiguity initialization can be significantly improved by robust schemes. Several robust schemes can effectively resist the simulated small and large gross errors.
Thirdly, to consider the positioning performance and computational efficiency, the IG-GIII scheme that includes the three-segment function is more effective than other proposed robust schemes in complex blocked environments.
Author Contributions: Methodology, software, and writing-original draft preparation, D.C.; conceptualization, validation, and writing-review and editing, Y.N.; investigation and formal analysis, S.W., W.S. and J.X.; investigation and data curation, J.B. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The data presented in this study are available on request from the authors on reasonable request.