Abstract

Precise point positioning (PPP) is used in many fields. However, pseudorange multipath delay is an important error that restricts its accuracy. Pseudorange multipath delay can be considered as the combination of effective information and observation noise; it can be modeled after removing the observation noise. In this work, elastic nets (EN) regularization denoising method is proposed and compared with L2-norm regularization denoising method. Then, quadratic polynomial (QP) model plus autoregressive (AR) model (QP + AR) are used to model the denoised pseudorange multipath delays. Finally, the modeling results are corrected to the observations to verify the improvement of BDS-3 single-frequency PPP accuracy. Three single-frequency PPP schemes are designed to verify the effectiveness of denoising method and QP + AR model. The experimental results show that the accuracy improvement of B3I and B2a is more obvious than that of B1I and B1C when the modeling values are corrected to the pseudorange observations. The improvement of B3I and B2a in the east (E) and up (U) directions can reach 10.6%∼34.4% and 5.9%∼65.7%, and the improvement of the north (N) direction is mostly less than 10.0%. The accuracy of B1I and B1C in E and U directions can be improved by 0%∼30.7% and 0.4%∼28.6%, respectively, while the accuracy of N direction can be improved slightly or even decreased. Using EN regularization denoising and QP + AR model correction, single-frequency PPP performs better at B3I and B2a, while L2-norm regularization denoising and QP + AR model correction perform better at B1I and B1C. The accuracy improvement of B2a and B3I is more obvious than that of B1I and B1C. The convergence time after MP correction of each frequency is slightly shorter. Overall, the proposed pseudorange multipath delays processing strategy is beneficial in improving the single-frequency PPP of BDS-3 satellite.

1. Introduction

BeiDou navigation satellite system (BDS) is a global satellite navigation system independently constructed and operated by China [1]. According to the “three-step” development strategy, the global service has been successfully launched in July 2020 [24]. As of June 2021, there are more than 50 BDS satellites in orbit, including BDS-2 and BDS-3 satellites, which are composed of GEO, MEO, and IGSO, respectively. Currently, five kinds of signals are transmitted by BDS satellites: B1I, B3I, B1C, B2a, and B2b. B1I and B3I are transmitted jointly by BDS-2 and BDS-3 satellites, which can ensure the smooth transition from BDS-2 to BDS-3 satellite [57]. B1C and B2a are transmitted by BDS-3 satellite, and the central frequency points are the same as GPS L1/L5 and Galileo E1/E5a, which can ensure better compatibility of BDS with GPS and Galileo [8]. In addition to providing BDS-3 navigation, positioning, and timing services, B2b can also provide services for international search and rescue and global short message communication [9, 10].

In the process of PPP using GNSS, many errors need to be corrected, such as antenna phase center correction and antenna phase wind-up. These errors can be corrected by the existing algorithms [1113]. However, some errors, such as pseudorange multipath delays, have not yet been corrected by an agreed upon method [14, 15]. The downlink navigation signals of satellites are transmitted under nonideal conditions, and the surrounding environment has a great influence on the received signals [2]. When the signal is reflected from the ground and surrounding objects into the receiver interior, the obtained observations contain pseudorange multipath delays and carrier phase multipath delays. The multipath delays of carrier phase observations are generally less than 1/4 wavelength, while the multipath delays of pseudorange observations can reach meter-level. Therefore, in some precision measurements, such as mining subsidence and building (structural) health monitoring, the correction of pseudorange multipath delays must be considered [15, 16]. The generation mechanism of pseudorange multipath delays is complex, and it has a large relationship with the surrounding reflection environment, so it is difficult to be accurately described [17]. However, the research shows that the pseudorange multipath delay is not a pure random error, but systematic to a certain extent. The systematic part can consider building models to correct the observations, which also makes it possible to deal with the pseudorange multipath delays.

At present, the modeling research of pseudorange multipath delays is mainly in time domain and space domain, and the models obtained by them take time and satellite position as independent variables, respectively. Time domain modeling is considering the periodicity of pseudorange multipath delays, so the application of sidereal filtering in pseudorange multipath delays is mainly studied [16, 18, 19]. Modeling in space domain is mainly considering the correlation between pseudorange multipath delays and satellite elevation, so the model construction of pseudorange multipath delays based on satellite elevation and azimuth has been studied [2, 20]. At present, BDS consists of BDS-2 and BDS-3, including GEO, MEO, and IGSO satellites with different orbital periods, which leads to difficulties in time domain modeling. In contrast, space domain modeling does not require so much periodicity of data and is more flexible. In addition, when the receiver is installed in some dynamic body, the reflection environment of the receiver is consistent, and the space domain method can also be used for modeling.

In 2012, Hauschild et al. found the pseudorange multipath delays in BDS-2 satellite signals for the first time [21]. Subsequently, relevant scholars studied the pseudorange multipath delays of BDS satellite but mostly studied the modeling of BDS-2 satellite, and the pseudorange multipath delays of BDS-3 satellite were less studied [2224]. It is generally believed that the pseudorange multipath delays of BDS-3 satellite are less than those of BDS-2 satellite, which has little influence on PPP. Therefore, the impact of pseudorange multipath delays of BDS-3 on PPP is rarely analyzed. By analyzing the pseudorange multipath delays of BDS-3, it is found that the pseudorange multipath delays of BDS-3 can reach decimeter-level or even meter-level under different epochs, which has a certain degree of influence on PPP [25, 26].

This paper mainly studied the pseudorange multipath delay of BDS-3 satellite, analyzed its characteristics, and modeled it based on the space domain, and the influence of correcting the pseudorange multipath delays of BDS-3 satellite on single-frequency PPP is analyzed. Firstly, the pseudorange multipath delays of BDS-2 and BDS-3 satellites were compared, and their characteristics were analyzed. Then, the EN regularization denoising method was proposed and compared with L2-norm regularization denoising. EN regularization adopts both L2-norm regularization and L1-norm regularization to maintain the original characteristics of the sequence but also has a certain sparsity [2730]. The pseudorange multipath delays after denoising were modeled by QP + AR model, and the modeling effects were analyzed. Finally, the modeling errors were corrected to the observations, and the positioning accuracy of BDS-3 satellite at each frequency was calculated.

The structure of the paper is arranged as follows. Section 2 presents the regularization denoising methods and the modeling principle of QP + AR model. Section 3 discusses the pseudorange multipath delays characteristics of BDS satellite, the regularization denoising effect, the modeling accuracy of QP + AR model, and the PPP accuracy of single-frequency. Finally, some conclusions and next step works are given in Section 4.

2. Methods

2.1. BDS Pseudorange Multipath Delays

In general, the linear combination of the original pseudorange observations and the carrier phase observations is used to construct multipath observations (MP) to evaluate the multipath delays characteristics of pseudorange observations [2, 5]. The MP of the s satellite can be expressed aswhere the subscripts i and j (i, j = 1, 2, 3; i ≠ j) represent different frequencies, P and represent pseudorange observations and carrier phase observations, f, λ denote frequency and wavelength, B is the combination of ambiguity term and hardware delay bias, N is the ambiguity term, d is a constant, and ε denotes other errors.

This combination of equations (1) and (2) can eliminate ionospheric delay, tropospheric delay, and geometric errors. When there is no cycle slip, B is generally a constant, which can be obtained by taking the average value of MP. Otherwise, segmented processing is required. When the constant term B is removed, the combined observations can be considered to contain only multipath delays and observation noise. In addition, the multipath delays of the carrier phase observation value is at most 1/4 cycle [16], and the impact is all in centimeter-level, which is negligible compared with the pseudorange multipath delays in decimeter-level or meter-level.

2.2. Regularization Denoising Method

The MP obtained in Section 2.1 contains pseudorange multipath delays and observation noise, so the observation noise should be removed before modeling. Common denoising methods include wavelet threshold denoising [31], empirical mode decomposition [32], and regularization denoising [15, 16]. In this work, EN regularization denoising method was proposed and compared with L2-norm regularization denoising. References [15, 16] confirmed that L2-norm regularization performed better in multipath delays denoising.

2.2.1. L2-Norm Regularization Denoising Method

As can be seen from the above content, MP contains pseudorange multipath delays and measurement noise, which can be specifically expressed aswhere represents the MP calculated in Section 2.1, represents the pseudorange multipath delays, represents the measurement noise, and n is the total number of epochs.

To extract pseudorange multipath delays signal, the target function can be constructed as follows:where ω represents the weight determined based on the elevation , , σc is the pseudorange accuracy, Ee is the satellite elevation of e epoch, μ is the Lagrange multiplier, and its solution method can be referred to in [15, 16, 33].

For the convenience of expression, equation (4) can be expressed in matrix form as follows:with

To obtain the optimal valuation of m, equation (5) should satisfy the following equation:

Hence, the MP sequence after L2-norm regularization denoising can be obtained as follows:

2.2.2. EN Regularization Denoising Method

According to Occam’s Razor principle, simpler models generally have better generalization performance [29]. Therefore, the concept of sparsity in machine learning is introduced here, and the sparse model will be simpler. When the L2-norm regularization is used, the solutions will tend to zero, but not equal to zero, while L1-norm regularization has higher penalty intensity, which can make some solutions equal to zero directly. Considering that L1-norm as a regular term can make the model have better sparsity [33], L1-norm is introduced with L2-norm to form EN regular term together. EN regularization combines the advantages of the two methods, but the computation is also relatively increased. Moreover, both L1-norm regular term and L2-norm regular term can compress the model coefficients. Therefore, EN regularization compresses the model coefficients twice, resulting in a large estimation deviation. However, Zou proposed that a scale transformation can be performed on the EN results to solve this problem [27].

Referring to equation (4), the target function of EN regularization can be constructed as follows:where μ1 and μ2 are the hyperparameters of L1-norm regular term and L2-norm regular term, respectively, and other parameters are the same as those in equation (4).

For the convenience of expression, equation (9) can be expressed in matrix form as follows:

According to equation (10), when μ1 = 1, μ2 = 0, the target function contains only L1-norm penalty term, and the EN regularization will be changed to L1-norm regularization. When μ1 = 0, μ2 = 1, the target function contains only L2-norm penalty term, and the EN regularization will be changed to L2-norm regularization.

Because the gradient of L1-norm is discontinuous, equation (10) cannot be solved directly [33]; it can be transformed as follows:where τ is a decimal slightly greater than zero.

Then, equation (10) can be expressed as follows:with

To obtain the optimal valuation of m, equation (12) should satisfy the following equation:

Hence, the MP sequence after EN regularization denoising can be obtained as follows:

However, due to the large dimension of the matrix in equation (10), the calculation time is increased. Through analysis, it can be seen that the coefficient matrix is an n-order tridiagonal matrix, which can be expressed as follows:with

In equation (15), the matrix elements satisfy the following conditions:

Therefore, equation (14) can also be solved by the Thomas algorithm with less computation [15].

2.3. MP Model Construction

The MP of BDS satellite has a certain correlation with the elevation, and the MP sequences after noise reduction can be modeled by QP model. The characteristics of each satellite are not completely consistent, so each satellite is modeled separately. The QP model of the i frequency of satellite s can be expressed as follows:where is calculated from equation (15), E(e) is the elevation of satellite s in the e epoch, ε is the model noise, and ai0, ai1, ai2 are the QP model coefficients, which can be solved by the weighted least square method.

Assuming that the observations of i frequency of S satellites are received in a certain period time, and the polynomial coefficients of all satellites are solved simultaneously, the matrix form of the error equation can be obtained as follows:with

After modeling with QP model, the residual presents a stable and random characteristic. AR model is used to further model the residual. The mathematical expression of AR model is as follows:where is the residual value of the satellite s of frequency i after QP modeling in epoch e, p is the order of AR model, is the AR coefficient of satellite s on the i frequency, and ε is the white noise.

In a certain period of time, can be solved by constructing the error observation equation as follows:with

The fitting accuracy of AR model is greatly affected by the order of the model. Akaike information criterion (AIC) and Bayesian information criterion (BIC) are commonly used to determine the order of AR model [34, 35]. Both AIC and BIC introduced penalty terms related to the number of model parameters to reduce the model parameters and avoid overfitting. However, AIC usually fails when the sample size is large, because the likelihood function value is large, which exceeds the influence of model parameters. In BIC, the punishment related to the number of model parameters is increased, which can effectively control the complexity of the model when the sample size is large. BIC is used to determine the model order, which can be expressed as follows [35, 36]:where is the mean square error of the model with , p is the number of model parameters, and n is the sample size. BIC value determines the quality of model, and the order with the smallest BIC value is the best order. Due to the large sample size and the amount of computation, the search range can be limited when searching for the optimal order. In this work, we set the order from 1 to 50 to find the optimal order.

3. Experiments and Analyses

In this section, the data processing method is verified by experiments. The observations of 17 iGMAS stations in 10 days (DOY 1-10, 2021) were collected, and all stations can receive BDS-3 and BDS-2 satellite signals. MEAN, RMS, and Range are selected as statistics to measure accuracy. Since there are few observations of B2b frequency, this work does not analyze its MP, nor study its single-frequency PPP. Figure 1 shows the distribution of iGMAS stations used in the experiment.

3.1. BDS MP Analysis

In this part, the MP of BDS satellites is mainly analyzed. WUH1 station (DOY 1) is selected as an example to visually display the MP characteristics of BDS-2 and BDS-3 satellites at different frequencies.

The MP characteristics of each satellite at each frequency are analyzed. In this section, one satellite of each type is selected as an example, which is shown in Figure 2. The missing epoch indicates the moment when the satellite cannot be observed. The values of blue, green, red, and cyan in the figure represent the statistical results of B1I, B3I, B1C, and B2a, respectively, and the black line represents the elevation change of the satellite. BDS-3 GEO satellites do not transmit B1C and B2a signals, so only MP at B1I and B3I frequencies is analyzed.

The experimental results show that each frequency of BDS satellite has MP, and it is inversely proportional to the satellite elevation. MEO and IGSO satellite have a large range of changes of elevation, while GEO satellite has a small range of changes, mostly within 5°. The MP of BDS-2 satellite changes with the elevation more obviously than that of BDS-3 satellite. The MP of different frequencies of the same satellite has obvious difference. For BDS-2 satellite, the influence of B3I is significantly less than that of B1I. Range and RMS show that the influence of B3I is more stable, and the fluctuation amplitude is smaller than that of B1I. For BDS-3 satellite, the fluctuation amplitude of MP at four frequencies shows the rule of B1C > B1I > B3I > B2a. The MP of B2a is relatively minimum, but it can still reach 0.5 m∼1.0 m at low elevation. B1C has the largest MP fluctuation amplitude. The MP of GEO satellite of BDS-2 and BDS-3 and MEO satellite of BDS-2 shows obvious periodicity at B1I and B3I frequencies. However, MEO and IGSO of BDS-3 did not show obvious periodicity. From the perspective of RMS and Range, the fluctuation amplitude of GEO satellite of BDS-2satellite is lower than that of MEO and IGSO satellite, while that of BDS-3 is opposite.

3.2. MP Denoising and Modeling

The regularization algorithm in Section 2.2 is used to denoise the MP sequence. Take C28 and C40 in Figure 2 as examples; both satellites contain four frequencies. For the convenience of distinction, the original MP sequence is denoted as MP, the L2-norm regularization denoising sequence is denoted as MP-regL2, and the EN regularization denoising sequence is denoted as MP-regEN. Figure 3 shows the comparison of the two satellites before and after the L2-regularization and EN regularization. Figure 4 shows the RMS and Range variations of the two satellites after MP regularization at each frequency. Table 1 shows the RMS statistical results of different satellites at different frequencies after regularization denoising.

It is obvious that the two regularization methods used in the experiment have the ability to denoise. The regularized sequence is still consistent with the volatility of the original MP sequence, but the sequence is relatively more stable. EN regularization makes the sequence sparse to a certain extent, and some small data in the original sequence become zero after EN regularization. In addition, the size of the hyperparameter is closely related to the stability of the sequence. The larger the hyperparameter is, the more stable the sequence is, but the data are prone to distortion. Therefore, it is extremely important to determine the size of the hyperparameter reasonably. Table 1 shows that the stability of MP sequences at different frequencies of different satellites is improved after regularization, and the stability after EN regularization is better than that after L2-norm regularization.

In this section, the modeling method in Section 2.3 is used to model MP. Most of the existing results suggest that the MP of BDS-2 MEO and IGSO satellites is significantly correlated with the elevation, while the MP of BDS-3 satellite is weakly correlated with the elevation, which makes it impossible to establish an effective model. In this work, QP + AR model is used to model the MP of BDS-3 satellite and analyze the modeling effect. Due to the small number of GEO satellites and the small range of elevation variation, only MEO and IGSO satellites of BDS-3 are studied in this part. The MP sequences after L2-norm regularization and EN regularization are modeled, respectively, and the model coefficients of each satellite at each frequency are calculated, respectively. Figures 5 and 6 show the modeling results for C28 and C40 in Figure 2.

Overall, the results of QP + AR model are basically consistent with the original sequence. The results of QP model mainly reflect the change of trend terms. In most cases, the impact on BDS-3 is within ±0.2 m, and the contribution to the modeling results is relatively small. By analyzing the coefficients of QP model, it is found that QP fitting parameters are relatively small. Table 2 lists the QP model parameters of C28 and C40. It is found that the coefficients of the QP model are small; especially, the coefficients of a2 are close to zero. AR model is more suitable for the modeling of stationary random sequences, which can accurately represent the fluctuation of sequences. However, it should be noted that the order of AR model is fixed during the modeling process. The order of AR adaptively selected by different epochs should be studied in the following work.

3.3. Single-Frequency PPP Experiment

In order to further test the modeling effect of Section 3.2, two iGMAS stations (WUH1 and KUN1; DOY 1-10, 2021) are selected for BDS-3 single-frequency PPP experiment. The precision ephemeris and clock bias products are mainly downloaded from MGEX. Table 3 shows some data processing strategies.

Three single-frequency PPP schemes are designed. Scheme 1: correct the errors without MP, and then carry out single-frequency PPP according to Table 3. Scheme 2: similar to scheme 1, but also correct MP. Using L2-norm regularization to denoise MP, the QP + AR is established, and the modeling results are corrected to the observation value, and then the single-frequency PPP is carried out according to Table 3. Scheme 3: similar to scheme 2, but EN regularization is used instead of L2-norm regularization.

To intuitively show the positioning accuracy of BDS-3 single-frequency PPP of the three schemes, the positioning results of WUH1 and KUN1 stations on DOY 1 are presented in Figures 7 and 8. Tables 4 and 5 show the positioning accuracy statistics of the two stations after convergence.

It can be seen that when BDS-3 satellite is used only for single-frequency PPP, the four frequencies in the E and N directions can reach centimeter-level after convergence, but the U direction is basically in the decimeter-level. When the observation environment of the station is consistent, the B2a frequency of BDS-3 performs best in single-frequency PPP.

Compared with Scheme 1, Scheme 2 and Scheme 3 show that the overall positioning accuracy is improved. Increasing MP correction can significantly improve the positioning accuracy of B3I and B2a frequencies, in which the improvement of the E and U directions can reach 10.6%–34.4% and 5.9%–65.7%, and the improvement of the N direction accuracy is mostly less than 10.0%. The positioning accuracy of B1I and B1C frequencies can be improved by 0%∼30.7% and 0.4%∼28.6% in the E and U directions, respectively, while the improvement in the N direction is small, and even the accuracy will be reduced.

Compared with Scheme 2, Scheme 3 performs better in the single-frequency PPP of B3I and B2a, while Scheme 2 performs better in the single-frequency PPP of B1I and B1C. The denoising of EN regularization is greater than that of L2-norm regularization, and it can make the sequence show a certain degree of sparsity. Therefore, the fluctuation amplitude of MP after denoising is lower than that of L2-norm regularization, which also leads to different QP + AR modeling values of the same epoch. As can be seen from Section 3.1, the MP of B1I and B1C is greater than that of B3I and B2a, so the MP correction of B1I and B1C frequency in Scheme 2 is greater than that in Scheme 3. The MP of B3I and B2a is relatively small, and the fluctuation is more stable, and in this case, the correction effect of Scheme 3 is better.

Convergence time is an important indicator of PPP performance. Generally, the positioning deviation in E, N, and U directions is less than 0.1 m and maintained for a certain time as the convergence condition [39]. In this work, BDS-3 single-frequency observation data is used for PPP, but the positioning accuracy is slightly lower, especially in the U direction. Based on the positioning results of many tests and considering the large positioning deviation in the U direction, the convergence condition in this work is defined as: the positioning deviation in the E and N directions is better than 0.1 m, in the U direction is better than 0.6 m, and the positioning deviation in the next 20 epochs is within the limit. Figure 9 shows the convergence time of WUH1 and KUN1 stations in DOY 1.

As can be seen from Figure 9, the convergence time of each frequency is quite different, among which B1I frequency has the longest convergence time, B3I frequency has the second, and B2a frequency has the shortest convergence time. In terms of convergence time and accuracy, the accuracy of the new signals transmitted by BDS-3 is better than that of the old signals, and the quality of B2a frequency observation data is relatively better. For the three schemes, the convergence time after MP correction is shortened to a certain extent, but it is not obvious. However, compared with Scheme 1 and Scheme 3, the convergence time of B1I frequency of KUN1 station is significantly shorter in Scheme 2. It is found that the positioning deviation of Scheme 1 and Scheme 3 has also changed smoothly after 7h, but the U direction is slightly greater than 0.6, so it is not counted as the convergence time. The convergence time of B1I frequency of WUH1 station is longer, which is also because the positioning deviation of the U direction is not less than 0.6 m. However, the positioning deviation of the three schemes also tends to be stable in about 6h. This is also the case for the positioning deviation of other frequencies. Due to the low accuracy of single-frequency PPP, how to set the convergence threshold rationally needs further study.

4. Conclusions

Pseudorange multipath delay is one of the main errors that affect the PPP, and an accurate correction model can improve the positioning accuracy. This work studied the influence of BDS-3 satellite MP on single-frequency PPP. Firstly, the characteristics of each frequency MP of BDS satellite are analyzed, and then the MP sequence is denoised by L2-norm regularization and EN regularization, respectively. The QP + AR model is constructed for the denoised sequence, and the modeling value is corrected to BDS-3 single-frequency PPP.

The following conclusions can be drawn from the experiments:(1)The MP of BDS-3 satellite can reach the decimeter-level. Although it is generally smaller than that of BDS-2 satellite, the influence cannot be ignored. Among the four frequencies, the MP of B1C is the largest, and that of B2a is the smallest. However, when the elevation is relatively low, the MP of BDS-3 can still reach 0.5 m∼1.0 m, which also indicates that the MP of BDS-3 is correlated with the elevation to a certain extent.(2)EN regularization is a combination of L1-norm regularization and L2-norm regularization. Compared with the L2-norm regularization denoising, MP sequences have certain sparsity and higher stability. However, it should be noted that the size of the regularization parameter has a great influence on the denoising effect, and its size needs to be determined reasonably.(3)The QP + AR is constructed by using the denoised MP and corrected to the pseudorange observations, which can improve the accuracy of single-frequency PPP. The accuracy improvement of B2a and B3I is more obvious than that of B1I and B1C, and the improvement in the E and U directions is significantly greater than that in the N direction. The accuracy of the E and U directions can be improved by 10.6%∼34.4% and 5.9%∼65.7%, and the N direction can be improved by less than 10.0%, but in some cases, the accuracy can even be reduced.(4)When B3I and B2a use EN regularization denoising sequence to construct QP + AR model to correct MP, their single-frequency PPP accuracy is better than that of L2-norm regularization denoising modeling to correct MP, while B1I and B1C use L2-norm regularization denoising modeling to correct MP with higher single-frequency PPP accuracy.(5)When BDS-3 single-frequency PPP converges, the position deviation in E and N directions can reach centimeter-level, and the U direction remains in demerit-level. The signal quality of B2a and B1C is better than that of B1I and B3I in terms of convergence time and accuracy.

In general, it is feasible to denoise MP and correct it by modeling, and the accuracy of BDS-3 single-frequency PPP can be improved to a certain extent by correcting MP. However, the number of stations used in this work is not large, which cannot be universal as the model parameters in reference [20]. In addition, the order of AR model is fixed when calculating model coefficients. In the subsequent research, it can be considered to change the order adaptively in different epochs. Although more computation is required, better modeling accuracy will be achieved. When the order of AR model is too large, we can also consider adding regularization to prevent overfitting. In addition, the positioning accuracy of single-frequency PPP is slightly lower, especially in the U direction, so how to reasonably determine the convergence threshold also needs to be further studied.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare no conflicts of interest.

Authors’ Contributions

Zhongchen Guo and Xuexiang Yu conceived the research scheme. Zhongchen Guo and Chao Hu verified the feasibility of the method and designed experimental scheme. Zhongchen Guo and Zhihao Yu implemented the software algorithms. Zhihao Yu and Chang Jiang collected the iGMAS data and analyzed data. Zhongchen Guo and Chao Hu wrote the manuscript.

Acknowledgments

The authors appreciate the International GNSS Monitoring and Assessment Service (iGMAS) and Multi-GNSS Experiment (MGEX) for the provision of relevant data and products. Thanks are due to Dr. Zhou Feng of Shandong University of Science and Technology for the open source software GAMP. This work was supported by the Open Research Fund of Coal Industry Engineering Research Center of Collaborative Monitoring of Mining Area’s Environment and Disasters (Grant No: KSXTJC202003), National Natural Science Foundation of China (Grant No: 41474026), Key Research and Development Program of Anhui Province (Grant No: 202104a07020014), Major Science and Technology Projects of Anhui Province (Grant No: 202103a05020026), Anhui Natural Science Foundation (Grant No: 2008085MD114), and Natural Science Foundation of Anhui Colleges (Grant no. KJ2020A0310).