Mitigation of Unmodeled Error to Improve the Accuracy of Multi-GNSS PPP for Crustal Deformation Monitoring

: High-rate multi-constellation global navigation satellite system (GNSS) precise point positioning (PPP) has been recognized as an e ﬃ cient and reliable technique for large earthquake monitoring. However, the displacements derived from PPP are often overwhelmed by the centimeter-level noise, therefore they are usually unable to detect slight deformations which could provide new ﬁndings for geophysics. In this paper, Global Positioning System (GPS), GLObalnaya NAvigatsionnaya Sputnikovaya Sistema (GLONASS), and BeiDou navigation satellite system (BDS) data collected during the 2017 Mw 6.5 Jiuzhaigou earthquake were used to further exploit the capability of BDS-only and multi-GNSS PPP in deformation monitoring by applying sidereal ﬁltering (SF) in the observation domain. The equation that uniﬁes the residuals for the uncombined and undi ﬀ erenced (UCUD) PPP solution on di ﬀ erent frequencies was derived, which could greatly reduce the complexity of data processing. An unanticipated long-term periodic error term of up to ± 3 cm was found in the phase residuals associated with BDS satellites in geostationary Earth orbit (GEO), which is not due to multipath originated from the ground but is in fact satellite dependent. The period of this error is mainly longer than 2000 s and cannot be alleviated by using multi-GNSS. Compared with solutions without sidereal ﬁltering, the application of the SF approach dramatically improves the positioning precision with respect to the weekly averaged positioning solution, by 75.2%, 42.8%, and 56.7% to 2.00, 2.23, and 5.58 cm in the case of BDS-only PPP in the east, north, and up components, respectively, and 71.2%, 27.7%, and 37.9% to 1.25, 0.81, and 3.79 cm in the case of GPS / GLONASS / BDS combined PPP, respectively. The GPS / GLONASS / BDS combined solutions augmented by the SF successfully suppress the GNSS noise, which contributes to the detection of the true seismic signal and is beneﬁcial to the pre- and post-seismic signal analysis. repeat for


Introduction
High-rate global navigation satellite systems (GNSSs) have shown great potential in observing both ground static and dynamic motions, which is more than a favorable complement to traditional seismometers [1,2]. Retrieving high-precision coseismic displacement in real time contributes The rest of the paper is organized as follows. The description of the multi-GNSS PPP data processing strategy is given first in Section 2. Afterward, the equations that describe the relationship between phase residuals on different frequencies are derived. Section 3 presents the multi-GNSS data description. Section 4 displays the analysis of the characteristics of BDS geostationary Earth orbit (GEO) phase residuals, PPP solutions using sidereally filtered multi-GNSS data, especially BDS PPP solutions, and a case study on the Jiuzhaigou earthquake. Finally, in Section 5, the conclusions and perspectives are provided.

PPP Model and Data Processing Strategy
After correcting the satellite orbits and clocks using the precise products, the PPP model for GPS, GLONASS, and BDS can be written as follows: r,j = −e G r · r r + t r + γ j,G I G r,1 + T G r + M G r,j + ξ G r,j p Rk r,j = −e Rk r · r r + t r + δ r,Rk + γ R j,Rk I Rk r,1 + T R r + M Rk r,j + ξ Rk r,j p C r,j = −e C r · r r + t r + δ r,C + γ j,C I C r,1 where p s r,j and l s r,j are the "observed minus computed" pseudo-range and phase observations from receiver r to satellite s on frequency j ( j = 1, 2); the superscripts G, R, and C represent the GPS, GLONASS, and BDS systems, respectively; R k denotes the frequency factor of the GLONASS satellite; e s r is the line-of-sight unit vector from receiver to satellite; r r denotes the vector of the receiver position increments with respect to the a priori position used for linearization; t r is the receiver clock offset; δ r,C and δ r,R k denote the inter-system bias (ISB) with respect to GPS and the inter-frequency bias (IFB) for GLONASS, respectively; and N s r,j refers to the float ambiguity containing code hardware and phase delay, while λ s j is its corresponding wavelength. The ionospheric delay at different frequencies can be expressed as I s r,j = γ j,s I s r,1 , γ j,s = f 2 1 / f 2 j ; T s r denotes the slant tropospheric delay; M s r,j and m s r,j represent the multipaths in pseudo-range and carrier-phase observations; and ξ s r,j and ε s r,j are the measurement noises. Other error items, such as the phase wind-up, relativity effects, Earth rotation, and tidal loading, are corrected by applying models described in Kouba [27]. It should be noted that, for GPS and GLONASS, the phase center offsets and variations (PCO and PCV) at both satellite and receiver are obtained from the International GNSS Service (IGS) antenna file. On the other hand, for BDS, the PCO and PCV at satellite are available from the IGS antenna file, but replaced by GPS at the receiver [28]. The data processing information is listed in Table 1 in detail. The receiver clock offset is treated as white noise, and estimated epoch by epoch. Together with other parameters, the receiver positions are estimated as daily solution using the least-squares estimator. Table 1. Data processing information for multi-global navigation satellite system (GNSS) precise point positioning (PPP).

Item Processing Information
Estimator Least-squares estimator for generating phase residuals Observations Raw pseudo-range and carrier-phase observations from GPS, GLONASS, and BDS Sampling rate 1 s Elevation cutoff 7 •

Weighting scheme
Elevation-dependent weight; 3 dm and 3 mm for GPS pseudo-range and carrier-phase; 4.5 dm and 3 mm for GLONASS pseudo-range and carrier-phase; 9.0 dm for BDS pseudo-range; and 5 mm and 15 mm for IGSO/MEO and GEO carrier-phase, respectively Satellite orbit/clock GBM final precise orbit/clock products generated by GFZ (Deng et al. [29]) Tropospheric delay The zenith hydrostatic delay corrected by Saastamoinen's model [30]; the zenith wet delay and the horizontal gradients estimated as piecewise constants every hour and six hours, respectively; For the SF method, the satellite orbital repeat times for GPS and BDS satellites are computed individually using Keplerian orbital elements from broadcast ephemerides [19,31]. On the other hand, for GLONASS, its broadcast ephemeris is presented by positions and velocities, and thus we used the aspect repeat time instead [32,33]. The phase residuals over n days before the day of interest are shifted by n times the orbital repeat time and then stacked for each station-satellite pair. Afterward, these stacked residuals are low-pass filtered with a cutoff frequency of 10 s to generate multipath corrections since any prominent multipath over shorter periods are not anticipated [34]. Finally, the phase observations on the day of interest are corrected, and then processed by the forward Kalman filter for simulating the kinematic PPP in real time.

Mathematic Relationship of Residuals on Different Frequencies
Since the coordinates are fixed and the zenith tropospheric delays and ambiguities are treated as constants during a period of time, the residuals, which are primarily the multipath errors, are related to the time-varying parameters, that is, receiver clock and ionospheric delays. Therefore, Equation (2) can be rearranged as follows: where A denotes the design matrix; E denotes a column vector of n-dimension with value one; I is an identity matrix of n-dimension; X is the parameter vector of receiver clock and ionospheric delays; L refers to the vector of unmodeled errors, and each element of l s j can be expressed as l s j + e s r · r r − T s r − λ s j N s j ; and V represents the vector of residuals. According to the least-squares criterion, the estimated parameters read as follows: By substituting Equation (5) into (3), the residuals are derived: The term A T A can also be simplified as follows: where a 11 = 2n, a 12 = a T 21 = −(1 + γ 2 )E T , and a 22 = (1 + γ 2 j,s )I. Applying the Gauss elimination method [35], the inverse matrix of A T A can be derived as follows: where Substituting Equation (8) into (6), the residual vector can be explicitly expressed as follows: Hence, the residuals of phase observations on different frequencies can be formulated as follows: which means the residuals on one frequency assimilate those on the other frequencies. Finally, we note that the relationship of residuals on two frequencies can be expressed as follows: This finding is of great significance, since we only need to calculate the multipath corrections on one frequency and can directly recover the corrections on another frequency by Equation (12). Note that in this paper, the residuals associated with one frequency are actually linearly combined residuals on different frequencies. A destructive Mw 6.5 earthquake occurred at 13:19:46 (UTC) on 8 August (DOY 220) 2017 in the Jiuzhaigou tourist area in Sichuan province of China at a relatively shallow depth of 20 km (http://news.ceic.ac.cn/). As shown in Figure 1, 14 GNSS stations from the Crustal Movement Observation Network of China (CMONOC) and BeiDou Ground Based Augmentation Systems (BDGBAS) networks were distributed near the epicenter. All the stations were capable of capturing GPS, GLONASS, and BDS signals with 1 Hz sampling rate. The data spans a period of time from DOY 206 to DOY 220. Generally, the orbital repeat times are 86,155 s for GPS, seven days and 84,442 s for GLONASS, 86,165 s for BDS GEO/IGSO, and six days and 84,697 s for BDS MEO, respectively. Taking the orbital repeat time and data length into consideration, the compromised number of days for residual stacking is seven for GPS and BDS IGSO/GEO, and one for GLONASS and BDS MEO.

Results and Discussion
The correctness of Equation (12) is validated first, followed by a detailed analysis of the GEO residuals. Then, the performance of BDS-only and multi-GNSS PPP is assessed. Finally, a case study of the Mw 6.5 Jiuzhaigou earthquake is shown.

Equation Validation
Figure 2 depicts the linear correlation between the phase residuals on two frequencies for GPS, GLONASS, and BDS at station SCPW on DOY 219. As can be seen, the residuals manifest strong negative correlation. All the correlation coefficients are -1.0, and the slopes of lines are -1.647, -1.652, and -1.673 for GPS, GLONASS, and BDS IGSO and MEO satellites, respectively, which are very close to the corresponding negative squares of ratios of the two frequencies (-1.648, -1.653, -1.672). Nonetheless, for BDS GEO satellites, the correlation coefficient is only about -0.4, and the slope of -0.458 shows a pronounced discrepancy with respect to the theoretical value of -1.672. This indicates that the phase residuals of the BDS GEO satellites on two frequencies have weak correlation and could not be properly described by Equation (12). The cause for this phenomenon is the pseudo-range bias, which degrades the precision of ionospheric parameters, consequently contaminating the residuals.

Results and Discussion
The correctness of Equation (12) is validated first, followed by a detailed analysis of the GEO residuals. Then, the performance of BDS-only and multi-GNSS PPP is assessed. Finally, a case study of the Mw 6.5 Jiuzhaigou earthquake is shown.

Equation Validation
Figure 2 depicts the linear correlation between the phase residuals on two frequencies for GPS, GLONASS, and BDS at station SCPW on DOY 219. As can be seen, the residuals manifest strong negative correlation. All the correlation coefficients are -1.0, and the slopes of lines are -1.647, -1.652, and -1.673 for GPS, GLONASS, and BDS IGSO and MEO satellites, respectively, which are very close to the corresponding negative squares of ratios of the two frequencies (-1.648, -1.653, -1.672). Nonetheless, for BDS GEO satellites, the correlation coefficient is only about -0.4, and the slope of -0.458 shows a pronounced discrepancy with respect to the theoretical value of -1.672. This indicates that the phase residuals of the BDS GEO satellites on two frequencies have weak correlation and could not be properly described by Equation (12). The cause for this phenomenon is the pseudo-range bias, which degrades the precision of ionospheric parameters, consequently contaminating the residuals. To investigate the pseudo-range weighting effects on carrier-phase residuals, the a priori precision of pseudo-range was set lower by a factor of two and three to 1.8 and 2.7 m, respectively. The results are shown in Figure 3. A higher a priori precision results in a higher observation weight. It is clear that the correlations rapidly increase from -0.4 to -0.9 with the decrease of pseudo-range weight. As expected, the corresponding slopes of the regression line become closer and closer to the square of the ratio of frequencies B1 and B2. Since the pseudo-range primarily provides the initial value for the least-squares estimator, to avoid contaminating the phase residuals, it is advisable to lower the weight of the pseudo-ranges. In this paper, the a priori precision of the BDS GEO pseudorange is set to 2.7 m.  To investigate the pseudo-range weighting effects on carrier-phase residuals, the a priori precision of pseudo-range was set lower by a factor of two and three to 1.8 and 2.7 m, respectively. The results are shown in Figure 3. A higher a priori precision results in a higher observation weight. It is clear that the correlations rapidly increase from -0.4 to -0.9 with the decrease of pseudo-range weight. As expected, the corresponding slopes of the regression line become closer and closer to the square of the ratio of frequencies B1 and B2. Since the pseudo-range primarily provides the initial value for the least-squares estimator, to avoid contaminating the phase residuals, it is advisable to lower the weight of the pseudo-ranges. In this paper, the a priori precision of the BDS GEO pseudo-range is set to 2.7 m. To investigate the pseudo-range weighting effects on carrier-phase residuals, the a priori precision of pseudo-range was set lower by a factor of two and three to 1.8 and 2.7 m, respectively. The results are shown in Figure 3. A higher a priori precision results in a higher observation weight. It is clear that the correlations rapidly increase from -0.4 to -0.9 with the decrease of pseudo-range weight. As expected, the corresponding slopes of the regression line become closer and closer to the square of the ratio of frequencies B1 and B2. Since the pseudo-range primarily provides the initial value for the least-squares estimator, to avoid contaminating the phase residuals, it is advisable to lower the weight of the pseudo-ranges. In this paper, the a priori precision of the BDS GEO pseudorange is set to 2.7 m.  To validate the correctness of Equation (12), the residuals on frequency two are recovered from frequency one, and then differenced with the observed residuals. The percentage error is based on the equation of (1.0 -a2/a1) × 100 %, where a1 and a2 are, respectively, the observed and projected residuals. As shown in Figure 4, the mean percentage errors for GPS, GLONASS, and BDS IGSO/MEO are all under 0.3 %, whereas for GEO they increase by about 6%-10 % to 2.5-4 mm. Although the recovery for GEO residuals performs not as well as that of other satellites, it is still acceptable, which gives a powerful proof of the correctness of Equation (12) To validate the correctness of Equation (12), the residuals on frequency two are recovered from frequency one, and then differenced with the observed residuals. The percentage error is based on the equation of (1.0 -a2/a1) × 100 %, where a1 and a2 are, respectively, the observed and projected residuals. As shown in Figure 4, the mean percentage errors for GPS, GLONASS, and BDS IGSO/MEO are all under 0.3 %, whereas for GEO they increase by about 6%-10 % to 2.5-4 mm. Although the recovery for GEO residuals performs not as well as that of other satellites, it is still acceptable, which gives a powerful proof of the correctness of Equation (12).

GEO Residual Analysis
Since the GEO satellites (C01-C05) are basically stationary relative to a point on the Earth's surface, the change rate of elevation angle is nearly zero, which means that the multipath from the ground should be close to a constant bias in the carrier-phase observation, theoretically. This is shown by Geng et al. [22], where the phase residuals of C01 satellite for about three hours at station CHPS are almost a constant with only a few fluctuations. To further investigate the characteristics of GEO residuals, the residuals at 27 globally distributed stations were calculated for 81 days. The station distribution map is depicted in Figure S1 (Supplementary Materials), and the information of five stations presented in this paper is listed in Table 2. Figure 5 typically delineates the residuals of C01 and C02 for about seven days. In contrast, excluding the daily peaks which are caused by the discontinuity of orbit and clock products at adjacent days, pronounced periodic errors of up to ± 3 cm can be found at some stations with the period of about a sidereal day. The residuals differ among

GEO Residual Analysis
Since the GEO satellites (C01-C05) are basically stationary relative to a point on the Earth's surface, the change rate of elevation angle is nearly zero, which means that the multipath from the ground should be close to a constant bias in the carrier-phase observation, theoretically. This is shown by Geng et al. [22], where the phase residuals of C01 satellite for about three hours at station CHPS are almost a constant with only a few fluctuations. To further investigate the characteristics of GEO residuals, the residuals at 27 globally distributed stations were calculated for 81 days. The station distribution map is depicted in Figure S1 (Supplementary Materials), and the information of five stations presented in this paper is listed in Table 2. Figure 5 typically delineates the residuals of C01 and C02 for about seven days. In contrast, excluding the daily peaks which are caused by the discontinuity of orbit and clock products at adjacent days, pronounced periodic errors of up to ± 3 cm can be found at some stations with the period of about a sidereal day. The residuals differ among stations for the same satellite. For example, the residuals of C01 and C02 for station GMSD contain subtle or even no periodic signal, which is consistent with the results provided by Geng et al. [22], whereas those for station CIBG reveal a conspicuous period. The residuals also differ among satellites for the same station. For example, for station JFNG, the residuals of C02 are totally different from those of C01 both in the amplitude and phase. This interesting periodic bias can be eliminated if differencing the observations between two nearby stations. That could be the reason why there are no relevant reports in the previous studies [23,36]. stations for the same satellite. For example, the residuals of C01 and C02 for station GMSD contain subtle or even no periodic signal, which is consistent with the results provided by Geng et al. [22], whereas those for station CIBG reveal a conspicuous period. The residuals also differ among satellites for the same station. For example, for station JFNG, the residuals of C02 are totally different from those of C01 both in the amplitude and phase. This interesting periodic bias can be eliminated if differencing the observations between two nearby stations. That could be the reason why there are no relevant reports in the previous studies [23,36].   Figure 6 presents the residual series of C01 and C03 from DOY 223 to 303, 2018 at stations DAE2, DAEJ, and GMSD, respectively. It is noteworthy that the amplitudes of residuals change rapidly on DOY 265 for C01 and on DOY 264 for C03, as the green dashed line shows. The residuals for C01 from DOY 223 to 264 and for C03 from DOY 263 to 303 are several times smaller than those on other days. Similar observations can also be made for C02, C04, and C05, as shown in Figure S2. Since these stations are hundreds of kilometers away from each other, it should not be associated with the environment. A preliminary conclusion is that this periodic bias originates from the satellite. The identification of the kind of bias is out of the scope of this paper, and will be investigated in the future.  Figure 6 presents the residual series of C01 and C03 from DOY 223 to 303, 2018 at stations DAE2, DAEJ, and GMSD, respectively. It is noteworthy that the amplitudes of residuals change rapidly on DOY 265 for C01 and on DOY 264 for C03, as the green dashed line shows. The residuals for C01 from DOY 223 to 264 and for C03 from DOY 263 to 303 are several times smaller than those on other days. Similar observations can also be made for C02, C04, and C05, as shown in Figure S2. Since these stations are hundreds of kilometers away from each other, it should not be associated with the environment. A preliminary conclusion is that this periodic bias originates from the satellite. The identification of the kind of bias is out of the scope of this paper, and will be investigated in the future.

Assessment of BDS-Only PPP with Multipath Correction
The a priori precision of GEO phase observations is usually set lower than GPS and GLONASS to account for the lower precision of orbit and clock products. However, the constant bias at the satellite can be partially absorbed by the time-invariant parameters such as ambiguities, and the remaining periodic errors can be mitigated by the SF approach. Therefore, it is expected that the positioning precision improves by properly setting the weight of the GEO phase observation. For this paper, four weighting schemes with a priori precisions of 15 mm, 10 mm, 6 mm, and 3 mm, named as (1) to (4), respectively, were designed. The smaller the precision value, the higher the observation weight, which contributes more to the PPP solution. For each scheme, two types of BDS-only solutions were calculated. For one solution, all the BDS observations were sidereally filtered. whereas for the other one, only the observations from the MEO and IGSO BDS satellites were filtered. About eight BDS satellites were visible on average during the observation, of which five were GEO satellites. Figure 7a shows the results without sidereal filtering for the GEO satellites at station GSWX on DOY 219, 2018, along with the mean root-mean-square (RMS) statistics for all stations. The weekly averaged positioning solutions are used as references. As can be seen, large wiggles, varying from a few centimeters to tens of centimeters, occur in all three components, especially for the up component. The best accuracies of the estimated displacements are achieved by weighting scheme (1) with the RMS values of 7.51, 3.50, and 11.38 cm for the east, north, and up components, respectively, and then they decrease with the increase of the a priori phase precision. The statistics are consistent with that of Li et al. [27]. This occurs because improperly setting the a priori precision higher magnifies the negative effects of the periodic errors in the GEO phases on the PPP solution.

Assessment of BDS-Only PPP with Multipath Correction
The a priori precision of GEO phase observations is usually set lower than GPS and GLONASS to account for the lower precision of orbit and clock products. However, the constant bias at the satellite can be partially absorbed by the time-invariant parameters such as ambiguities, and the remaining periodic errors can be mitigated by the SF approach. Therefore, it is expected that the positioning precision improves by properly setting the weight of the GEO phase observation. For this paper, four weighting schemes with a priori precisions of 15 mm, 10 mm, 6 mm, and 3 mm, named as (1) to (4), respectively, were designed. The smaller the precision value, the higher the observation weight, which contributes more to the PPP solution. For each scheme, two types of BDS-only solutions were calculated. For one solution, all the BDS observations were sidereally filtered. whereas for the other one, only the observations from the MEO and IGSO BDS satellites were filtered. About eight BDS satellites were visible on average during the observation, of which five were GEO satellites. Figure 7a shows the results without sidereal filtering for the GEO satellites at station GSWX on DOY 219, 2018, along with the mean root-mean-square (RMS) statistics for all stations. The weekly averaged positioning solutions are used as references. As can be seen, large wiggles, varying from a few centimeters to tens of centimeters, occur in all three components, especially for the up component. The best accuracies of the estimated displacements are achieved by weighting scheme (1) with the RMS values of 7.51, 3.50, and 11.38 cm for the east, north, and up components, respectively, and then they decrease with the increase of the a priori phase precision. The statistics are consistent with that of Li et al. [27]. This occurs because improperly setting the a priori precision higher magnifies the negative effects of the periodic errors in the GEO phases on the PPP solution. After applying sidereal Remote Sens. 2019, 11, 2232 11 of 20 filtering on GEO carrier-phase observations, the displacement noises are effectively alleviated for all schemes, as shown in Figure 7b.
As the optimal results among four schemes, scheme (2) has the smallest RMS values, which are dramatically reduced by 75.2%, 42.8%, and 56.7%, compared to the unfiltered ones, to 2.00, 2.23, and 5.58 cm for the three components, respectively. This indicates that the SF approach can significantly improve the precision of GEO phase observations to around 10 mm, thereby making a better contribution to the PPP solution. Consequently, the a priori phase precision of GEO satellites was set to 10 mm for the following experiments. After applying sidereal filtering on GEO carrier-phase observations, the displacement noises are effectively alleviated for all schemes, as shown in Figure 7b. As the optimal results among four schemes, scheme (2) has the smallest RMS values, which are dramatically reduced by 75.2%, 42.8%, and 56.7%, compared to the unfiltered ones, to 2.00, 2.23, and 5.58 cm for the three components, respectively. This indicates that the SF approach can significantly improve the precision of GEO phase observations to around 10 mm, thereby making a better contribution to the PPP solution. Consequently, the a priori phase precision of GEO satellites was set to 10 mm for the following experiments. The power spectral densities (PSDs) for scheme (2) were also calculated using Welch's method for each station, and then averaged for all the PSDs for the specific frequency from all stations, as shown in Figure 8. It is clear that the SF approach mainly reduces the PSDs on the longer periods over 2000 s. In other words, the periodic errors of the GEO phase observations primarily weigh on the lowest frequency band. The PSD reductions are 3 dB and 5 dB, for the north and up components. In contrast, for the east component, the reduction can reach up to 20 dB since the periodic errors are primarily projected to the east-west direction due to the special distribution of the GEO satellites. The power spectral densities (PSDs) for scheme (2) were also calculated using Welch's method for each station, and then averaged for all the PSDs for the specific frequency from all stations, as shown in Figure 8. It is clear that the SF approach mainly reduces the PSDs on the longer periods over 2000 s. In other words, the periodic errors of the GEO phase observations primarily weigh on the lowest frequency band. The PSD reductions are 3 dB and 5 dB, for the north and up components. In contrast, for the east component, the reduction can reach up to 20 dB since the periodic errors are primarily projected to the east-west direction due to the special distribution of the GEO satellites.

Assessment of Multi-GNSS PPP with Multipath Correction
As revealed by Geng et al. [16], the integration of GPS and GLONASS data led to a substantial reduction of the high-rate displacement noise by up to 40% compared to a GPS-only solution within European regions. In this section, the two major purposes investigated are presented. One shows the benefits from multi-GNSS in improving the precision of displacements; the other one concerns the question of whether the multi-GNSS data can suppress the effect of GEO periodic error. Four experiments numbered from (1) to (4) were designed, and the corresponding results of station GSWX in the east, north, and up direction are displayed in Figures 9a-d, 10a-d, and 11a-d, respectively. Table 3 lists the data filtering strategy for the four experiments. Experiment (1) gives the results without sidereal filtering, whereas Experiments (2) to (4) present the sidereally filtered results. Note that the GEO observations are not filtered in Experiment (2), and they are excluded in Experiment (4).  Figure 9a, it is observed that although the combined GPS/GLONASS/BDS (GRC) solution shows about 12.7% improvement compared to the GC combination, both of them are evidently inferior to the GPS-only solution. The fusion of GRC can effectively reduce most noise in contrast to the BDS solution without sidereal filtering, leading to a decline in the RMS values from 7.51, 3.50, and 11.38 cm to 4.34, 1.12, and 6.10 cm in the three components, respectively (see Figure 10; Figure  11), but there are still many fluctuations over the period of 5000 s remaining. After the sidereal filtering, a considerable reduction of 21.9% and 23.5% in terms of RMS can be found in the east and up components (see Figure 11), respectively, for GPS, whereas the SF approach may occasionally introduce undesirable ramps as shown in the north component around 6 to 7 h (see Figure 10). From a comparison between Figures 9a and 9b, the application of SF on BDS IGSO and MEO observations can slightly reduce them by about 3.2 dB noise over periods from 50 s to 5000 s, but it fails in alleviating noise over longer periods. The RMS values for sidereally filtered GC and GRC solutions are comparable to those without sidereal filtering, which are 4.91 and 4.13 cm, respectively. Once the

Assessment of Multi-GNSS PPP with Multipath Correction
As revealed by Geng et al. [16], the integration of GPS and GLONASS data led to a substantial reduction of the high-rate displacement noise by up to 40% compared to a GPS-only solution within European regions. In this section, the two major purposes investigated are presented. One shows the benefits from multi-GNSS in improving the precision of displacements; the other one concerns the question of whether the multi-GNSS data can suppress the effect of GEO periodic error. Four experiments numbered from (1) to (4) were designed, and the corresponding results of station GSWX in the east, north, and up direction are displayed in Figure 9a-d, Figure 10a-d, and Figure 11a-d, respectively. Table 3 lists the data filtering strategy for the four experiments. Experiment (1) gives the results without sidereal filtering, whereas Experiments (2) to (4) present the sidereally filtered results. Note that the GEO observations are not filtered in Experiment (2), and they are excluded in Experiment (4). From Figure 9a, it is observed that although the combined GPS/GLONASS/BDS (GRC) solution shows about 12.7% improvement compared to the GC combination, both of them are evidently inferior to the GPS-only solution. The fusion of GRC can effectively reduce most noise in contrast to the BDS solution without sidereal filtering, leading to a decline in the RMS values from 7.51, 3.50, and 11.38 cm to 4.34, 1.12, and 6.10 cm in the three components, respectively (see Figure 10; Figure 11), but there are still many fluctuations over the period of 5000 s remaining. After the sidereal filtering, a considerable reduction of 21.9% and 23.5% in terms of RMS can be found in the east and up components (see Figure 11), respectively, for GPS, whereas the SF approach may occasionally introduce undesirable ramps as shown in the north component around 6 to 7 h (see Figure 10). From a comparison between Figure 9a and 9b, the application of SF on BDS IGSO and MEO observations can slightly reduce them by about 3.2 dB noise over periods from 50 s to 5000 s, but it fails in alleviating noise over longer periods. The RMS values for sidereally filtered GC and GRC solutions are comparable to those without sidereal filtering, which are 4.91 and 4.13 cm, respectively. Once the GEO observations are sidereally filtered (Experiment (3)), as delineated in Figure 9c, Figure 10c, and Figure 11c, the RMS values of the GRC solutions are reduced to 1.25, 0.81, and 3.79 cm from 4.34, 1.12, and 6.10 cm of Experiment (1), for the east, north, and up components, respectively, which are dramatic improvements of about 71.2%, 27.7%, and 37.9%, respectively. The PSDs across almost the entire frequency band decline substantially, especially for the longer periods where the PSDs decline on average by about 10.2 dB. Excluding the GEO satellites can improve the precision of GC and GRC solutions to some extent; however, it also increases the RMS values by 14.8% and 13.6% to 1.63 and 1.42 cm, respectively, compared with the sidereally filtered counterparts. Moreover, it seriously deteriorates the precision of the BDS-only PPP because of the poor geometry of the satellite constellation. Additionally, the sidereally filtered BDS-only solution (se Figure 7b) outperforms the unfiltered GPS solution in the east and up components (see Figures 9a and 11a), but is slightly worse than the sidereally filtered GPS solution in the north and up components (see Figures 10b and 11b).  Figure 7b) outperforms the unfiltered GPS solution in the east and up components (see Figures 9a and 11a), but is slightly worse than the sidereally filtered GPS solution in the north and up components (see Figures 10b and 11b).

A Case Study for the Mw6.3 Jiuzhaigou Earthquake
Two stations, SCJZ and GSZQ, close to the epicenter were selected to carry out the experiments, and the corresponding results are shown in Figure 12; Figure 13. The final coseismic displacements for these two stations were calculated using 15 days of GPS data before and 4 days after the event from [37]. For station SCJZ, although the data were interrupted after 48,032 s, the seismic waveforms were well recorded. The largest amplitudes of about 4 cm occurred in the north component for all the types of solutions. Compared with the north component for GPS, the seismic signal in the east displacement was indistinguishable from the noise, which might mislead preseismic analysis. In contrast, for BDS, the signal-to-noise ratio improved with the smaller fluctuations ahead of the event.

A Case Study for the Mw6.3 Jiuzhaigou Earthquake
Two stations, SCJZ and GSZQ, close to the epicenter were selected to carry out the experiments, and the corresponding results are shown in Figure 12; Figure 13. The final coseismic displacements for these two stations were calculated using 15 days of GPS data before and 4 days after the event from [37]. For station SCJZ, although the data were interrupted after 48,032 s, the seismic waveforms were well recorded. The largest amplitudes of about 4 cm occurred in the north component for all the types of solutions. Compared with the north component for GPS, the seismic signal in the east displacement was indistinguishable from the noise, which might mislead preseismic analysis. In contrast, for BDS, the signal-to-noise ratio improved with the smaller fluctuations ahead of the event. The GRC solution without sidereal filtering was biased by about −1.30 cm in the east component, which was reduced by 89.2% to -0.14 cm after filtering. For station GSZQ, the largest fluctuation happened in the east component with a peak value of about -2.3 cm. The seismic waveforms embedded in the vertical direction were seemingly overwhelmed by high-level noise. Affected by the remaining systematic errors, some fluctuations of 2 to 3 cm spanning over several minutes could still be found in the east component for the BDS solutions. Because of the data interruption at station SCJZ, 2 h displacements of the sidereally filtered GRC solution before and after the arrival time of seismic waves were used to estimate the static offsets of station GSZQ, which were 1.9 mm and 5.4 mm with respect to the references 0.4 ± 1.2 mm and 3.6 ± 0.8 mm, in the east and north components, respectively. Overall, the superiority of the GRC solution with sidereal filtering over a single-system or unfiltered solution in alienating low-frequency errors on tens of seconds to minutes is clearly demonstrated. The GRC solution without sidereal filtering was biased by about −1.30 cm in the east component, which was reduced by 89.2% to -0.14 cm after filtering. For station GSZQ, the largest fluctuation happened in the east component with a peak value of about -2.3 cm. The seismic waveforms embedded in the vertical direction were seemingly overwhelmed by high-level noise. Affected by the remaining systematic errors, some fluctuations of 2 to 3 cm spanning over several minutes could still be found in the east component for the BDS solutions. Because of the data interruption at station SCJZ, 2 h displacements of the sidereally filtered GRC solution before and after the arrival time of seismic waves were used to estimate the static offsets of station GSZQ, which were 1.9 mm and 5.4 mm with respect to the references 0.4 ± 1.2 mm and 3.6 ± 0.8 mm, in the east and north components, respectively. Overall, the superiority of the GRC solution with sidereal filtering over a single-system or unfiltered solution in alienating low-frequency errors on tens of seconds to minutes is clearly demonstrated.   The black, blue, and purple lines refer to the sidereally filtered G, C, and GRC solutions, respectively, while the green lines refer to the GRC solution without sidereal filtering.

Conclusions
For this study, the performance of multi-GNSS PPP, especially BDS PPP, in monitoring subtle deformation was investigated. The sidereal filtering approach was employed to mitigate the multipath in phase observations to improve the precision of PPP. The equations describing the relationship between phase residuals on different frequencies were rigorously derived, which could significantly reduce the complexity of multipath processing. A satellite-dependent periodic error term with an amplitude of up to ± 3 cm was found in the BDS GEO phase residuals, one of the main sources that limits the precision of BDS PPP. The results indicated that the systematic errors originated mainly from GEO, had periods longer than 2000 s, and could not be alleviated by the fusion of multi-GNSS, whereas the multipath errors from IGSO and MEO had periods from 50 to 5000 s. Traditionally, GEO observations are weighted lower to account for the imprecise orbit and clock products, but this also reduces the GEO's contribution to the solution. The SF approach can effectively mitigate the periodic errors, thus improving the precision of the GEO phase to around 10 mm. Compared with the BDS-only PPP solutions without sidereal filtering, the one using the SF Figure 13. Displacements at station GSZQ during the Jiuzhaigou earthquake. The lines are shifted vertically to avoid overlap. The yellow dashed horizontal lines denote the mean of displacements. The black, blue, and purple lines refer to the sidereally filtered G, C, and GRC solutions, respectively, while the green lines refer to the GRC solution without sidereal filtering.

Conclusions
For this study, the performance of multi-GNSS PPP, especially BDS PPP, in monitoring subtle deformation was investigated. The sidereal filtering approach was employed to mitigate the multipath in phase observations to improve the precision of PPP. The equations describing the relationship between phase residuals on different frequencies were rigorously derived, which could significantly reduce the complexity of multipath processing. A satellite-dependent periodic error term with an amplitude of up to ± 3 cm was found in the BDS GEO phase residuals, one of the main sources that limits the precision of BDS PPP. The results indicated that the systematic errors originated mainly from GEO, had periods longer than 2000 s, and could not be alleviated by the fusion of multi-GNSS, whereas the multipath errors from IGSO and MEO had periods from 50 to 5000 s. Traditionally, GEO observations are weighted lower to account for the imprecise orbit and clock products, but this also reduces the GEO's contribution to the solution. The SF approach can effectively mitigate the periodic errors, thus improving the precision of the GEO phase to around 10 mm. Compared with the BDS-only PPP solutions without sidereal filtering, the one using the SF approach can effectively improve the positioning accuracy, with respect to the weekly averaged positioning solutions, by 75.2%, 42.8%, and 56.7% to 2.00, 2.23, and 5.58 cm in the east, north, and up components, respectively. It is comparable to that of GPS in the east component, and slightly worse in the north and up components. After applying sidereal filtering, the accuracy of the combined GPS, GLONASS, and BDS solution can also be improved by 71.2%, 27.7%, and 37.9% to 1.25, 0.81, and 3.79 cm in the three directions, respectively, compared to the unfiltered results.
Author Contributions: K.Z. and X.Z. conceived and designed the experiments; K.Z. and P.L. performed the research; K.Z. and X.Z. wrote the paper; X.L., X.C., J.S., M.G., and H.S. provided advice and reviewed the paper.