#### Abstract

Conventional ambiguity resolution (AR) strategy of directly fixing all raw frequency ambiguities for uncombined precise point positioning (PPP) generally requires a long initial time to ensure fixing reliability. The rich signals from multi-frequency GNSS can provide new opportunities for rapid PPP AR. In this paper, based on a unified multi-frequency uncombined PPP model, the raw frequency ambiguities are linearly transformed to the between-satellite-single-differenced extra-wide-lane (EWL), wide-lane (WL), and narrow-lane (NL) ambiguities, and a cascading ambiguity resolution (CAR) method of fixing EWL/WL/NL ambiguities sequentially is proposed. Meanwhile, a partial ambiguity fixing (PAF) strategy with ambiguity subset adaptively selected based on the successively increased elevations is also adopted in each step to improve the fixing rate. Further, experiments with globally distributed stations are carried out to verify the algorithm. With the constraints of EWL/WL AR, the precision of NL ambiguity and its variance-covariance matrix can be effectively optimized, so the time to first fix (TTFF) is significantly shortened, and the ratio value is also improved to varying degrees. As for BDS-only solution, the average TTFF is shortened from 23.7 min to 10.9 min, with an improvement over 50%; for GPS/BDS, GPS/Galileo, and GPS/Galileo/BDS combined solutions, the TTFF is shortened from 14.3 min, 9.4 min, and 6.7 min to 8.1 min, 2.6 min, and 1.8 min, which are, respectively, shortened by 43.4%, 72.3%, and 73.1%. In general, the proposed CAR strategy can shorten the TTFF of multi-GNSS multi-frequency PPP to about 2 min. The performance of EWL/WL/NL ambiguity-fixed solutions is also analyzed. The NL solution is generally at centimeter-level accuracy over the entire period; however, limited by the atmosphere errors during the convergence stage, the EWL/WL solutions can only obtain decimeter-level accuracy, and the difference between them and NL solution gradually decreases with the continuous improvement of atmosphere accuracy.

#### 1. Introduction

Precise point positioning (PPP) technology can achieve unified high-precision positioning on a global scale, and it does not require the support of dense reference stations, which greatly improves the flexibility of Global Navigation Satellite System (GNSS) positioning and thus has been widely used in geodetic and geodynamic applications [1–3]. However, due to the weak model strength and strong correlation among various parameters, centimeter-level accuracy often requires a convergence time of about 10∼20 min [4, 5]. Ambiguity resolution (AR) can shorten the convergence time to a certain extent and improve the positioning reliability at the same time [6]. The core of successful AR is to eliminate the fractional cycle bias (FCB) in float ambiguity and then restore its integer characteristics. In recent years, a wide range of studies have been studied concerning reliable FCB estimation for integer ambiguity-fixing [7–9]. Objectively speaking, real-time PPP AR has now become feasible for precise positioning users. This is of great significance to promote the application of PPP in many fields. However, in the real-time perspective, it still needs significant improvement since it cannot provide almost instantaneous positioning likewise the widely used real-time kinematic (RTK) positioning technology, such as standard RTK (SRTK) and network RTK (NRTK).

For speeding up the PPP AR, a well-known approach is to provide external atmospheric corrections for each satellite through accurate regional atmospheric modelling [10, 11]. Then, the tropospheric and ionospheric errors can be corrected directly instead of being estimated as unknown parameters. This enhanced PPP recently attracts a lot of research attention, since it can achieve almost equivalent positioning performance as NRTK. However, in some special areas, e.g., the marine area, this mode is unrealistic as no dense reference network can be provided for inversing the atmospheric corrections. Another approach for speeding up the PPP positioning is using multi-GNSS data. With the continuous development and improvement of the Global Positioning System (GPS), the Global Navigation Satellite System (GLONASS), the Galileo and the BeiDou Navigation Satellite System (BDS), and the new generation of navigation satellites can all broadcast triple- or multi-frequency signals, and GNSS has officially entered an era of multi-system coexistence. The multi-system and multi-frequency also provide new opportunities for the improvement of positioning performance.

The benefit with the increased number of satellites can be interpreted as enhanced satellite geometry and the improved AR fixing rate with a large number of satellites. From another point of view, as many satellites can broadcast triple-frequency, quad-frequency, and even five-frequency signals, the initialization time of PPP is also expected to be shortened with multi-frequency observations. Based on simulated GPS triple-frequency data, Geng and Bock construct a triple-frequency PPP model, which uses ambiguity-fixed-ionosphere-free (AFIF) wide-lane (WL) observable instead of the raw pseudorange to tightly constrain the position parameters, and the initialization time of narrow-lane (NL) ambiguity is greatly shortened [12]. Compared with dual-frequency case, the convergence time can be slightly shortened by adding a few GPS Block IIF satellites [13, 14]. Using real-time BDS-2 triple-frequency data, Gu et al. assisted the fixing of B1 ambiguity through the resolved two WL ambiguity of B1-B2 and B2-B3, and the results show that the contribution of the third frequency is not obvious due to the poor BDS-2 geometry [15]. Li et al. also adopted a similar approach to analyze the AR performance of BDS-2 triple-frequency PPP in the Asia-Pacific region, and the results show that the performance of PPP in limited environments can be improved with the addition of the third frequency [16]. Li et al. combined BDS-2 B1I/B2I/B3I and Galileo E1/E5a/E5b triple-frequency to speed up the PPP performance by resolving the extra-wide-lane (EWL) and WL ambiguities and the NL ambiguities from ionosphere-free (IF) combination [17]. Geng et al. improved the triple-frequency PPP model by using more multi-frequency signals from GPS/BDS-2/Galileo/QZSS [18]. In their research, raw ambiguities are mapped at the normal equation level into their EWL, WL, and NL counterparts for integer cycle resolution. By implementing the single-epoch EWL and WL AR, Geng et al. also conduct the instantaneous decimeter-level positioning with the EWL/WL model [19, 20]. According to those studies, we can see that the multi-constellation and multi-frequency signals can effectively contribute to the faster convergence of PPP. It can be expected that, with more frequencies involved in AR and positioning calculation, better PPP performance can be obtained.

In this study, we exploit the fast PPP model with cascading ambiguity resolution, using the integrated multi-frequency data from GPS, Galileo, BDS-2, and BDS-3. In each AR step, the partial AR strategy based on successively increased elevations is adopted to improve the AR fixing rate. We will analyze the benefits in detail that the PPP AR speed and positioning can be improved using the multi-GNSS multi-frequency cascading fixing strategy.

#### 2. Materials and Methods

In this section, we describe the triple-frequency and quad-frequency observation equations and corresponding PPP AR model. The used multi-frequency signals of GPS, Galileo, BDS-2, and BDS-3 are shown in Table 1.

##### 2.1. Multi-Frequency PPP Observation Equations

The raw observation equations for undifferenced pseudorange and carrier phase can be expressed as follows:where indices refer to the satellite, receiver, and carrier frequency band, respectively; and are pseudorange and phase observations; is the geometric range between receiver and satellite; and are zenith wet-troposphere delay (ZWD) and the corresponding mapping function; is the speed of light in vacuum; and are receiver and satellite clock error, respectively; is line-of-sight ionosphere delay at the first frequency; is the ionosphere factor; is the integer ambiguity with wavelength ; and are receiver and the satellite code bias; and are receiver and satellite phase bias; and denote unmodelled errors in pseudorange and phase observations.

Based on the assumption that the satellite hardware bias is stable over time, it is usually lumped with ambiguity and estimated as constants, and this also makes the ambiguities nonintegers. With the development of multi-GNSS and multi-frequency, some scholars have found that, for some specific types of satellites, such as GPS Block IIF and BDS-2, there are obvious periodic changes in satellite phase bias, and then the concept of inter-frequency clock bias (IFCB) is introduced [21–23]. Meanwhile, the phase bias is further written as follows [24]:where and denote the time-invariant and time-varying parts of satellite phase bias, respectively.

Due to linear dependence between various estimated parameters, (1) is rank deficient. To get a full-ranked function model, the reparameterization process is usually carried out. Firstly, satellite-related dual-frequency IF combination of code bias and time-varying parts of phase bias is absorbed by precise clock, and similarly, the receiver-related dual-frequency IF combination of code bias is lumped with receiver clock bias:where can be expressed as linear combination of the hardware delays on raw frequency:with .

Secondly, in order to eliminate the correlation between the ionosphere and code bias, and at the same time to compensate for the time-varying parts of phase bias, the following reparameterization is also required:where DCB denotes differential code bias (DCB), and for a given frequency and , its definition is as follows:

Finally, with the above reparameterization process, the full-ranked observation model is obtained as follows:where and denote receiver inter-frequency bias (IFB) and satellite DCB; among them, the satellite DCB is usually corrected using products released by analysis centers in advance [25]; denotes float ambiguity; denotes the residual phase bias, and since the weight of pseudorange is low, it does not need to be considered additionally [24, 26]; and denotes IFCB. The specific expressions are as follows:

The estimated parameters of the above multi-frequency uncombined PPP model include

##### 2.2. Cascading PPP AR Model

With the above full-ranked multi-frequency uncombined PPP observation equations, we can obtain the ambiguities on each raw carrier-phase observation. After being applied with satellite FCB corrections and satellite-differenced (SD) operation, the produced SD ambiguities will have the integer property theoretically. Assume that the so-called float solution of the multi-frequency GNSS model together with its variance-covariance matrix (vc-matrix) is denoted aswhere the superscripts represent GPS, Galileo, BDS-2, and BDS-3, respectively; represents the SD ambiguities on each frequency; is the float solution of the parameter vector except ambiguities, including the coordinate components, the zenith tropospheric delay, and the slant ionospheric delay on each satellite. For the multi-frequency signals on GPS L1/L2/L5, Galileo E1/E5a/E6/E5b, BDS-2 B1I/B2I/B3I, and BDS-3 B1c/B1I/B2a/B3I, the ambiguity vector for each system can be expressed as

The wavelengths of the original ambiguities in (11) are ranging from 19.03 cm to 25.48 cm. Affected by observation noise and unmodelled errors, sometimes, these original ambiguities may be not able to be fixed reliably. Generally, the longer the wavelength of the ambiguity, the easier it is to be fixed, since the effect of unmodelled errors can be reduced. Therefore, we can transform those original ambiguities to the combination form with longer wavelengths, such as EWL or WL. Based on (10), the estimated unknown parameters and the corresponding vc-matrix can be transformed using (12) and (13):

In (12) and (13), the subscripts , , represent NL, WL, and EWL, respectively. By (12) and (13), the original ambiguities can be transformed to the ambiguities with different wavelengths. In related RTK positioning studies, sequential fixing of EWL, WL, and NL ambiguities will perform better as the previously fixed ambiguities increase the parameter precision for later fixings. Here, referring to the previous studies for the combination choices [27, 28], we adopt the combination classification in Table 2 as the EWL, WL, and NL selections. The corresponding fixing sequence of ambiguities is shown in Figure 1. It needs to be noted that there are two EWL combinations in quad-frequency cases, i.e., Galileo and BDS-3, and these two kinds of EWL ambiguities will be fixed synchronously, instead of sequentially.

After determining the EWL, WL, and NL combinations, the conversion matrix in (12) and (13) can be specifically expressed as

In the fixing process, the three kinds of ambiguities will be fixed in descending order of wavelengths. Firstly, the EWL ambiguities are fixed to their integers using the least-squares ambiguity decorrelation adjustment (LAMBDA) method [29]. The EWL ambiguities have long wavelengths and usually can be fixed easily with large ratio values, even with a single epoch. After the EWL ambiguities are fixed, the remaining parameters including WL ambiguities, NL ambiguities, and other nonambiguity parameters can be updated using the following equation:

And the corresponding vc-matrix can be similarly updated as

In (15) and (16), the elements with the superscript “” denote the EWL ambiguity-fixed parameter estimation. With the constraint of ambiguity-fixed EWL ambiguities, the precisions of float WL ambiguities will be improved. Then, sequentially, the constrained WL ambiguities are also fixed as with LAMBDA method. Similarly with (15) and (16), the remaining parameters including NL ambiguities and other nonambiguity parameters can be further updated with the fixed WL ambiguities, as (17) and (18)

In (17) and (18), the elements with the superscript “” denote the WL ambiguity-fixed parameter estimation. Based on (17) and (18), the new NL ambiguities are further fixed as . Then, the final parameters can be updated with (19),where the elements with the superscript “” denote the final NL ambiguity-fixed parameter estimation.

Compared with the conventional ambiguity resolution model for directly fixing the original ambiguities, the above described cascading AR (CAR) strategy can fully use beneficial constraints of fixed ambiguities from the previous step. Therefore, the AR performance, especially the time to first time (TTFF), can be improved theoretically. It needs to be noted that the prerequisite for the operability of the above CAR strategy is that the ambiguities in each step need to be fixed correctly. The research by Deo and El-Mowafy [30, 31] pointed out that the EWL/WL ambiguities can be reliably fixed using the geometry-free (GF) model. However, due to the large pseudorange noise, multiple epochs of data are usually required to smooth the noise so as to ensure the fix reliability. Thus, in this paper, we resolve all the EWL, WL, and NL ambiguities with the geometry-based (GB) model in each step. During the data processing, in order to ensure the AR reliability, the resolved optimal integer candidates at each step will be validated with ratio test with the threshold of 2.0, instead of the integer rounding (IR) strategy.

##### 2.3. Partial Ambiguity Fixing Strategy in Each AR Step

It is generally believed that the low-elevation satellites are easily affected by observation noise, multipath, and residual atmospheric errors, and their ambiguities accuracy is relatively poor, so it might be difficult to fix all ambiguities simultaneously [32]. Although the LAMBDA algorithm can give the optimal solution for all ambiguities, it may not pass the subsequent ambiguity acceptance test due to the influence of the low elevation ambiguities. In addition, for multi-frequency uncombined PPP with high-dimensional ambiguity, it is not necessary to fix all ambiguities when the number of satellites is sufficient. Therefore, the partial ambiguity fixing (PAF) strategy, with a subset of ambiguities fixed instead of all ambiguities, is proposed. The core of PAF lies in the selection of ambiguity subsets, and the commonly used indicators include elevation, precision of ambiguity, and continuous tracking number. In this paper, the ambiguity subset is selected based on the successively increased elevations [33]. In each step of CAR process, all ambiguities are first to be fixed; if the ratio test is passed, the corresponding ambiguity-fixed solution can be obtained; otherwise, the PAF is performed with specific process as follows: Step 1: sort the ambiguities of EWL/WL/NL in ascending order of elevation; Step 2: based on a predefined starting elevation (for example, 10°), filter the qualified ambiguity subsets, and then fix it through the LAMBDA algorithm. If the Ratio value meet the user-defined threshold, the acceptance test is passed; otherwise, proceed to Step 3 Step 3: increase the starting elevation in Step 2 by 5°, and then repeat Step 2. In each fixing process, if the number of available ambiguities in the subset is less than 5 or the starting elevation is larger than 50°, the PAF procedure will be terminated and keep float solution at current epoch.

In addition to the above models and strategies, other specific processing options for multi-frequency PPP are listed in Table 3.

#### 3. Experiments and Results

##### 3.1. Data Description

The multi-constellation and multi-frequency observation data of GPS, Galileo, and BDS from the International GNSS Service (IGS) Multi-GNSS Experiment (MGEX) network and Australian Regional GPS Network (ARGN), as well as the rapid precise satellite orbit and clock correction provided by GeoForschungsZentrum Potsdam (GFZ) are used to validate the PPP performance with the proposed method. It is worth mentioning that, due to the lack of precise products for some satellites, the number of available BDS-3 satellites is actually 23. Firstly, observation data on 2020 DOY 148 from 303 stations with 30 s sampling intervals are used for FCB estimation [7, 8], and the distribution of these stations is shown in Figure 2. These stations cover a variety of mainstream receiver types. Then, multi-frequency observation data from 18 stations with 30 s sampling intervals are used for positioning validation, and the distribution of these stations is shown in Figure 3. All the 18 stations can receive the GPS/BDS-2 triple-frequency and Galileo/BDS-3 quad-frequency signals. The detailed site information of these stations, including the site name, receiver type, antenna type, and average satellite number per epoch, is summarized in Table 4. During the process, the Geostationary Earth Orbit (GEO) satellites of BDS-2 and BDS-3 are excluded due to the poor accuracy of their observations and orbit products.

##### 3.2. Experiment of BDS (BDS-2 and BDS-3) Multi-Frequency Observation

Firstly, taking BDS-2 and BDS-3 multi-frequency data as example, we analyze the multi-frequency PPP performance. With the station GAMG as the case study, the number of visible BDS-2 and BDS-3 satellites and the corresponding Position Dilution of Precision (PDOP) are shown in Figure 4. In the period, the average visible satellite numbers of BDS-2 and BDS-3 are 9.5 and 5.2, respectively. We can see that, with BDS-2 and BDS-3, the available satellites and the observation structure are enough for the PPP process. The ratio values of EWL, WL, and NL ambiguity fixing with LAMBDA method are shown in Figure 5. We can see that the ratio values of EWL and WL are obviously larger than those of NL. This indicates that the EWL and WL ambiguity fixing are much easier than NL.

As in Figure 6, the time to first fix (TTFF) of the PPP processing at station GAMG is 17.5 min. To better illustrate, the subfigure shows the positioning errors for the first half hour. The criterion of the first ambiguity fixing is that the ratios reach the threshold, and the subsequent 10 epoch can also satisfy the Ratio condition. After the EWL, WL, and NL ambiguities are fixed, centimeter-level positioning can be obtained, as in Figure 6. The root mean square (RMS) of positioning accuracies in north direction (N), east direction (E), and up direction (U) is 1.0 cm, 1.0 cm, and 3.8 cm, respectively.

In order to compare the convergence time of the conventional uncombined AR (UC-AR) and the proposed CAR, the TTFFs of the six stations in Asia-Pacific area (within red dotted box in Figure 3) are shown in Figure 7(a). From Table 3, we can see in the Asia-Pacific area that the visible BDS satellites are more than other areas, especially for the BDS-2 satellites. The average TTFF with UC-AR of the six stations is 23.7 min, while, with CAR, the average TTFF is reduced to 10.9 min, with an improvement over 50%. All ratios of NL ambiguity fixing with the two methods are shown in Figure 7(b). We can see that the ratios with CAR method are obviously larger than those of UC-AR.

##### 3.3. Experiment of Multi-GNSS Integration

Besides the above BDS-2 and BDS-3 satellites, GPS and Galileo multi-frequency observations can be further combined to improve the PPP performance. With the station KOUR at South America as the case study, the numbers of visible satellites with GPS/BDS-2/BDS-3, GPS/Galileo, and GPS/Galileo/BDS-2/BDS-3 are shown in Figure 8(a). The GPS satellites with only dual-frequency signals were also used, which participated in positioning calculation from the second step, i.e., the WL AR. We can see that the combination of GPS/BDS-2/BDS-3 and the combination of GPS/Galileo have the comparable number of visible satellites, which ranges from about 15–20. The combinations of GPS/Galileo/BDS-2/BDS-3 have the most visible satellites, which ranges from about 20–30. The positioning performance with the above three combinations is shown in Figure 8(b), including the TTFF and the positioning accuracies. The TTFFs for the three combinations are 11.0 min, 2.5 min, and 2.0 min. The combination of GPS/BDS-2/BDS-3 consumes the longest TTFF, which may be caused by less visible satellites at the initialization phase. It is worth noting that the TTFF with GPS/Galileo/BDS-2/BDS-3 is only 2.0 min. The convergence speed is dramatically improved compared with the traditional dual-frequency PPP, which usually costs 10–30 min. The positioning accuracies with ambiguity-fixed solutions for the three combinations have no significant differences, since there are abundant visible satellites for each combination.

**(a)**

**(b)**

The statistics of the TTFFs of all the 18 stations in Table 3 are shown in Figure 9, where the performance of UC-AR and CAR is presented in contrast. The results of GPS/BDS-2/BDS-3, GPS/Galileo, and GPS/Galileo/BDS-2/BDS-3 are shown in the upper, middle, and lower panels, respectively. We can see that whether, for any combination of the three cases, the CAR method can obviously reduce the TTFFs, compared with the UC-AR method, particularly, with the combination of GPS/Galileo/BDS-2/BDS-3, almost all the stations except station BREW can realize the first fix within 2 min, and the average TTFF over the 18 stations is only 1.8 min. This means that, with the multi-constellation and multi-frequency observations, one can realize the rapid PPP solutions within 2 min all over the world. All ratios of NL ambiguity fixing with the two methods are shown in Figure 10. Similar to the conclusion in Figure 7, we can see that the ratios with CAR method are also obviously larger than those of UC-AR.

For a more detailed analysis for the UC-AR method and CAR method, we further compare the precision of NL ambiguities from some satellites at station KIRU. The ambiguity precision is extracted from the vc-matrix, which can reflect the model strength of AR model. The results of each satellite pair from each system are shown in Figure 11. It can be clearly seen that, compared UC-AR method, the precision of NL ambiguities is improved to varying degrees with constraints of EWL/WL AR previously, and its convergence speed is faster. To a certain extent, this also leads to higher search efficiency and shorter TTFFs for subsequent NL AR.

Different form the UC-AR method, the CAR method can simultaneously obtain EWL, WL, and NL ambiguity-fixed solutions, so the performance of different ambiguity-fixed solutions is further compared. Take the results of GODS station as an example, as shown in Figure 12, and the subfigures also show the positioning errors for the first 15 min. The TTFF of NL solution is 2 min, while the EWL and WL solutions can be instantaneously fixed. Benefiting from the low noise level of NL combination, its accuracy is generally at centimeter-level. In contrast, for EWL and WL solutions, since the accuracy of atmosphere parameters is low during the convergence stage, even if the ambiguities are fixed successfully, their accuracy is still poor, and with the continuous improvement of atmosphere accuracy in the filter process, the deviation between them and NL solution gradually decreases.

#### 4. Discussion

In uncombined PPP AR process, the conventional UC-AR strategy of directly fixing all raw frequency ambiguities is a theoretically optimal method; however, due to the short wavelength of NL ambiguity, it usually takes a certain convergence time to achieve the desired accuracy, so as to ensure the reliability of ambiguity fixing. In addition, compared with the traditional dual-frequency IF PPP model, the dimension of multi-frequency uncombined PPP is high, and it is extremely time-consuming to search and fix it as a whole. In contrast, in the proposed CAR strategy, the linearly transformed EWL, WL, and NL ambiguities are fixed and updated sequentially to reduce the dimension of fixed ambiguities in each step, which is expected to effectively improve the AR efficiency. At the same time, considering that the significance of ratio-test decreases with the increase of ambiguity dimension [39, 40], so for the same judgment criteria, the CAR strategy is expected to achieve higher ratio values and epoch fix rates. Although, in the CAR strategy, the NL ambiguity still needs multiple epoch of data to achieve reliable fixing, its precision and vc-matrix are optimized under the constrains of previous EWL and WL AR, thus achieving a shorter TTFF. We also found that when the atmosphere accuracy is sufficient, the difference between the EWL/WL/NL solutions is small, as shown in Figure 12, and this phenomenon inspires us to further utilize the atmosphere augmentation generated by regional reference stations and adopt WL ambiguity-fixed solution, instead of NL, to provide precise positioning with several centimeters. After all, the WL ambiguity is easier to be fixed, even in single-epoch mode [19, 20].

#### 5. Conclusions

In this paper, we focus on fast AR strategy of multi-GNSS and multi-frequency uncombined PPP. Based on a brief introduction of the multi-frequency PPP mathematical model, a CAR method of fixing linearly transformed EWL/WL/NL ambiguity sequentially is proposed, and at the same time, the specific process of PAF is also described. In order to verify the algorithm, experiments are carried out using globally distributed MGEX stations. In terms of AR performance, compared with the conventional UC-AR method, the CAR strategy can optimize the precision of NL ambiguity and its vc-matrix through EWL/WL AR, and the TTFFs are significantly shortened. As for BDS-only solution, the average TTFFs are shortened from 23.7 min to 10.9 min, with an improvement over 50%; for GPS/BDS, GPS/Galileo and GPS/Galileo/BDS combined solutions, the TTFFs are shortened from 14.3 min, 9.4 min, and 6.7 min to 8.1 min, 2.6 min, and 1.8 min, which are, respectively, shortened by 43.4%, 72.3%, and 73.1%. In general, the CAR strategy can shorten the TTFFs of multi-GNSS multi-frequency PPP to about 2 min. In addition, the CAR strategy can also increase the NL AR ratio values to varying degrees. Since the CAR strategy can obtain EWL/WL/NL solutions simultaneously, the performance of different types of ambiguity-fixed solutions is further analyzed. Among them, the NL solution has the highest accuracy, generally at centimeter-level over the entire period; however, the EWL and WL solutions are affected by atmosphere errors during the initial convergence stage, and the accuracy is about several decimeters. With the continuous improvement of atmosphere accuracy in the filter process, the deviation between them and NL solution gradually decreases.

#### Data Availability

The multi-GNSS precise product provided by GFZ is available at ftp://ftp.gfz-potsdam.de/pub/GNSS/products/mgex/ (accessed on 1 November 2021). The MGEX DCB product provided by Chinese Academy of Sciences (CAS) is available at https://cddis.nasa.gov/archive/gnss/products/bias/ (accessed on 1 November 2021). The MGEX observation data from IGS is available at https://cddis.nasa.gov/archive/gnss/data/daily/ (accessed on 1 November 2021). The MGEX observation data from ARGN is available at ftp://ftp.data.gnss.ga.gov.au/daily/ (accessed on 1 November 2021).

#### Conflicts of Interest

The authors declare no conflicts of interest.

#### Acknowledgments

The authors sincerely thank IGS, ARGN, GFZ, and CAS for providing multi-GNSS data and products. This research was funded by the National Natural Science Foundation of China (Grant nos. 41904022 and 41774027) and the Fundamental Research Funds for the Central Universities (2242021R41134).