Corrections of BDS Code Multipath Error in Geostationary Orbit Satellite and Their Application in Precise Data Processing

Multipath error is a main error source in Global Navigation Satellite System (GNSS) data processing, which cannot be removed by a differential technique because of the strong relationship with the environment around the station. The multipath effect of the code observables is more complex than that of the carrier-phase observables, especially for BeiDou Navigation Satellite System (BDS) geostationary orbit (GEO) satellites. In this contribution, we deeply analyzed the characteristic and effect on the precise data processing of GEO satellite multipath errors based on a large number of permanent GNSS stations. A linear combination of code and carrier-phase observables was used to analyze the characteristics of repeatability for BDS GEO’s multipath. Then, a correction method was proposed to eliminate the multipath error of the GEO code observables, based on wavelet transform. The experiment data were collected at 83 globally distributed stations, from multi-GNSS experiments and national BDS augmentation systems, from days 32 to 66 in 2017. The results show that the systematic multipath variation component of the GEO code observables can be obtained with wavelet transform, which can significantly contribute to correcting the multipath error of GEO satellites. The average root mean square error (RMSE) of the multipath series is decreased by approximately 19.5%, 20.2%, and 7.5% for B1, B2, and B3, respectively. In addition, some experiments, including ionospheric delay extraction and satellite clock estimation, were conducted in simulated real-time mode in order to validate the effect of the correction methods. For the ionospheric delay estimation, the average RMSE of the slant ionospheric delay is reduced by approximately 15.5%. Moreover, the multipath correction can contribute greatly to shortening the convergence time of the satellite clock estimation of the BDS GEO satellites.


Introduction
Multipath is a major error in Global Navigation Satellite System (GNSS) error sources. As multipath is associated with the site environment and cannot be cancelled out by differential technique, it will exert considerable effects on the quality of GNSS signals, and should be well-treated in high-precision GNSS data processing. For example, severe multipath effects will slow down the convergence of precise point positioning (PPP) and decrease the reliability of PPP ambiguity resolution [1,2]. As for its positioning in a dynamic environment, sometimes, a poor geometry results in the positioning accuracy being more sensitive to multipath errors [3][4][5]. As for the ionospheric delay estimation,

Mathematical Models
In this section, the basic mathematical models, including the multipath combination method and wavelet transform method, are briefly introduced.

Multipath Combination
Assuming that the code rang and the carrier phase observables of the receiver (r), to the visible satellite (s) at a certain frequency (f i ) are P s r,i and Φ s r,i , which can be expressed by the following equation: P s r,i = ρ s r + t r,sys − t s + b r,i − b s i + trop s r + iono s r,i + ε P Φ s r,i = ρ s r + t r,sys − t s + trop s r − iono s where ρ s r denotes the geometric distance from the station to the satellite; t r,sys denotes the receiver clock error associated with GNSS; t s denotes the satellite clock error; b r,i and b s i denote the hardware delay deviation of the code rang observables in frequency (f i ) for the receiver and satellite, respectively; trop s r and iono s r,i denote the tropospheric delay and ionospheric delay on the signal propagation path, respectively; N s r,i is the float ambiguity, consisting of integer ambiguity and initial phase bias; λ i is the phase wavelength; and ε P and ε Φ denote the noise on the code range and carrier phase observables, respectively. The multipath combination is a linear combination of single-frequency code range and dual-frequency carrier-phase observables [28], which can be described as follows: where C s r, f i is a multipath combination; i and j denote different frequencies of f i and f j , respectively; and i j; The multipath combination derived by Equation (2) contains not only the multipath, but also the ambiguity, hardware delay, and noise. When no cycle slip occurs, the ambiguity is assumed to be a constant value, and can be derived through averaging over the epochs. Thus, a zero-mean term for multipath effect on frequency f i can be derived as follows: where MP s r,i is a zero-mean term multipath combination, and C s r,i is the average value of C s r,i over time. Theoretically, multipath on carrier-phase observables will not exceed one quarter of the carrier wave cycle, which is approximately several centimeter's [18]. In comparison with the multipath on code observables, the multipath on carrier-phase observables can be neglected. Thus, multipath combination can be used for code multipath analysis. However, the multipath series derived through Equation (3) consists of low-frequency and high-frequency components, which should be separated in order to analyze the characteristic of multipath errors.

Wavelet Transform
Wavelet transform is used to represent or approach a signal with a family of wavelet functions generated from a prototype function. It has been widely used in GNSS data processing [14,21,29,30]. According to research, GPS bias terms, such as multipath and ionospheric delay, behave like low-frequency noise, and the observation noise behaves as high-frequency noise. Hence, the GPS bias terms are concentrated in the narrow low-frequency band, and a high-frequency resolution is needed in order to identify them [30]. Multi-resolution analysis can provide a formal approach to constructing the wavelet basis. The basic concept of multi-resolution analysis is to analyze the signal at different scales by using filters of different cut-off frequencies. The signal is passed through a series of high-pass filters (HPF) in order to analyze the high frequencies, and it is passed through a series of low-pass filters (LPF) in order to analyze the low frequencies. Thus, wavelet transform can be used to achieve enough frequency resolution to discriminate the multipath error in the original GNSS observables [13]. Figure 1 illustrates the multi-resolution analysis process using wavelet transform. Applying a narrow daughter wavelet to the original signal is equivalent to applying a HPF, which completes path 1. Extracting the leading low-frequency requires applying a number of daughter wavelets that are wider than the signal you need to match, then applying a final daughter wavelet that becomes a HPF completing path 2. Moreover, the low-frequency components obtained by the decomposition can be used for signal reconstruction, and detailed information can be found in the literature [21]. As for the details of the wavelet transform, they are listed in the Appendix A.
Sensors 2019, 19, x FOR PEER REVIEW 4 of 15 Figure 1 illustrates the multi-resolution analysis process using wavelet transform. Applying a narrow daughter wavelet to the original signal is equivalent to applying a HPF, which completes path 1. Extracting the leading low-frequency requires applying a number of daughter wavelets that are wider than the signal you need to match, then applying a final daughter wavelet that becomes a HPF completing path 2. Moreover, the low-frequency components obtained by the decomposition can be used for signal reconstruction, and detailed information can be found in the literature [21]. As for the details of the wavelet transform, they are listed in the Appendix. The key issue of wavelet transform is the determination of the wavelet function and decomposition level, but there is no unified theoretical standard so far. The Daubechies wavelets are a good option, as they are orthonormal bases of compactly supported wavelets [31]. We did not carry out any specific theoretical analysis on the selection of wavelet basis functions, but directly used the Daubechies wavelet on the basis of a previous study [32]. The decomposition level was determined to be three for data processing.

Experimental Analysis
Given that the multipath error is related to the station environment, if the station environment remains the same, then the difference in multipath error among continuous days will be extremely small. Therefore, the multipath error can be extracted from the historical data, and can be applied to the subsequent data processing. In this section, 83 globally distributed stations were used for the multipath characteristic analysis of BDS GEO satellites. Some of these stations were used for their ionospheric delay extraction and satellite clock estimation in simulated real-time mode in order to verify the effectiveness of the proposed method.

Data Collection and Processing Strategy
A GNSS network with 83 stations, which were from the multi-GNSS experiment campaign (MGEX) of the international GNSS service (IGS) and the National BDS Augmentation Service System (NBASS), was used in this paper [33][34][35]. Among them, the data of MGEX was achieved through the Crustal Dynamics Data Information System (CDDIS) [36]. The station distribution is shown in Figure  2. Some stations could not track the GEO satellites, which are mostly located in the area of America. Only the observables from 78 stations could be used for GEO satellite analysis. The data were collected from days 32 to 66 of 2017 with 30 s intervals, and the corresponding calendar dates were from 1 February 2017 to 7 March 2017.
The multipath combination consists of multipath errors and noises of observables, which can be divided into low-frequency and high-frequency components. These random noises have a low or no The key issue of wavelet transform is the determination of the wavelet function and decomposition level, but there is no unified theoretical standard so far. The Daubechies wavelets are a good option, as they are orthonormal bases of compactly supported wavelets [31]. We did not carry out any specific theoretical analysis on the selection of wavelet basis functions, but directly used the Daubechies wavelet on the basis of a previous study [32]. The decomposition level was determined to be three for data processing.

Experimental Analysis
Given that the multipath error is related to the station environment, if the station environment remains the same, then the difference in multipath error among continuous days will be extremely small. Therefore, the multipath error can be extracted from the historical data, and can be applied to the subsequent data processing. In this section, 83 globally distributed stations were used for the multipath characteristic analysis of BDS GEO satellites. Some of these stations were used for their ionospheric delay extraction and satellite clock estimation in simulated real-time mode in order to verify the effectiveness of the proposed method.

Data Collection and Processing Strategy
A GNSS network with 83 stations, which were from the multi-GNSS experiment campaign (MGEX) of the international GNSS service (IGS) and the National BDS Augmentation Service System (NBASS), was used in this paper [33][34][35]. Among them, the data of MGEX was achieved through the Crustal Dynamics Data Information System (CDDIS) [36]. The station distribution is shown in Figure 2. Some stations could not track the GEO satellites, which are mostly located in the area of America. Only the observables from 78 stations could be used for GEO satellite analysis. The data were collected from days 32 to 66 of 2017 with 30 s intervals, and the corresponding calendar dates were from 1 February 2017 to 7 March 2017. correlation between different days. Therefore, the two parts must be separated, and the multipath error must be extracted from the original multipath series.

GEO Satellite Multipath Calibration and Elimination
In this section, the characteristics of the GEO satellite multipath on triple-frequency observables are analyzed and discussed. Then, the low-frequency part of the multipath series is extracted and applied to correct multipath error.

Multipath Calibration
BDS GEO satellites can provide triple-frequency code observables. Thus, an analysis is performed to illustrate whether the multipath errors at different frequencies are correlated. Figure 3 presents the one-week multipath time series on a triple-frequency of BDS GEO satellite (C04) at station CUT0. To understand the characteristics of the multipath series more clearly, the multipath series at triple-frequency of date-of-year (DOY) 37 (6 February 2017) were also amplified separately. For simplicity, multipath combinations on the frequency B1, B2 and B3 are abbreviated to MP1, MP2 and MP3, respectively. As is shown in Figure 3, the amplitude fluctuations of the multipath series vary with the frequency. For example, the amplitude of MP1 series is about 2.0 m, which is significantly larger than those of MP1 and MP3. Overall, the multipath series of all of the triplefrequencies shows evident periodic repeatability. The repetition period is approximately one day, which is close to the GEO satellite orbit period. On the basis of the above analysis, the periodical features might be used to eliminate the multipath effects and to improve the quality of the precise data processing. The multipath combination consists of multipath errors and noises of observables, which can be divided into low-frequency and high-frequency components. These random noises have a low or no correlation between different days. Therefore, the two parts must be separated, and the multipath error must be extracted from the original multipath series.

GEO Satellite Multipath Calibration and Elimination
In this section, the characteristics of the GEO satellite multipath on triple-frequency observables are analyzed and discussed. Then, the low-frequency part of the multipath series is extracted and applied to correct multipath error.

Multipath Calibration
BDS GEO satellites can provide triple-frequency code observables. Thus, an analysis is performed to illustrate whether the multipath errors at different frequencies are correlated. Figure 3 presents the one-week multipath time series on a triple-frequency of BDS GEO satellite (C04) at station CUT0. To understand the characteristics of the multipath series more clearly, the multipath series at triple-frequency of date-of-year (DOY) 37 (6 February 2017) were also amplified separately. For simplicity, multipath combinations on the frequency B1, B2 and B3 are abbreviated to MP1, MP2 and MP3, respectively. As is shown in Figure 3, the amplitude fluctuations of the multipath series vary with the frequency. For example, the amplitude of MP1 series is about 2.0 m, which is significantly larger than those of MP1 and MP3. Overall, the multipath series of all of the triple-frequencies shows evident periodic repeatability. The repetition period is approximately one day, which is close to the GEO satellite orbit period. On the basis of the above analysis, the periodical features might be used to eliminate the multipath effects and to improve the quality of the precise data processing. The method of wavelet decomposition and reconstruction proposed is applied so as to derive the low-frequency component of the multipath series. With the wavelet decomposition and reconstruction, the low-frequency and high-frequency components in the multipath time series for each frequency can be separated. Figure 4 shows the three-level Daubechies wavelet decomposition of the multipath series. As can be seen from Figure 4, the extrema and amplitudes of the noise decrease with the increase of the decomposition level. Moreover, the low-frequency component extracted form the third level is extremely close to the original multipath series. So, it is also called the approximation. Thus, the low-frequency components are expected to be used for correcting the original observables of both the current and subsequent periods.  The method of wavelet decomposition and reconstruction proposed is applied so as to derive the low-frequency component of the multipath series. With the wavelet decomposition and reconstruction, the low-frequency and high-frequency components in the multipath time series for each frequency can be separated. Figure 4 shows the three-level Daubechies wavelet decomposition of the multipath series. As can be seen from Figure 4, the extrema and amplitudes of the noise decrease with the increase of the decomposition level. Moreover, the low-frequency component extracted form the third level is extremely close to the original multipath series. So, it is also called the approximation. Thus, the low-frequency components are expected to be used for correcting the original observables of both the current and subsequent periods.  The method of wavelet decomposition and reconstruction proposed is applied so as to derive the low-frequency component of the multipath series. With the wavelet decomposition and reconstruction, the low-frequency and high-frequency components in the multipath time series for each frequency can be separated. Figure 4 shows the three-level Daubechies wavelet decomposition of the multipath series. As can be seen from Figure 4, the extrema and amplitudes of the noise decrease with the increase of the decomposition level. Moreover, the low-frequency component extracted form the third level is extremely close to the original multipath series. So, it is also called the approximation. Thus, the low-frequency components are expected to be used for correcting the original observables of both the current and subsequent periods.    . "Raw" represents the original multipath series, "LF" represents the lowfrequency component extracted from the multipath series, and "Residual" indicates the multipath residual series after subtracting the low-frequency component.

Multipath Elimination
According to Section 3.2.1, the low-frequency component of the multipath series derived by the wavelet transform can greatly decrease the effect of the multipath error. To evaluate the improvement due to the low-frequency removal, the root mean square error (RMSE) of the one-week multipath series was computed for both without and with the multipath correction, respectively. Figure 6 shows the RMSE of the GEO satellite (C04) at station CUT0. The blue bars represent the RMSE without the multipath error correction, while the red bars represent the RMSE with the multipath error correction. It is obvious that the RMSE of the multipath series considerably decreases with the low-frequency components of the previous day, being removed from the raw time series. The RMSE of the multipath series at station CUT0 on the frequencies of B1, B2, and B3 is decreased by 42.0%, 34.0%, and 60.0% according to the average RMSE, respectively. The RMSE of MP3 is smaller than that of MP1 and MP2, mainly because of its low noise of observable. Table 1 shows the statistical results of the multipath series without and with the multipath correction for each GEO satellite from all of the stations. With the multipath correction, the accuracy of the multipath series at three frequencies is improved. For example, the RMSE of C01 satellite is decreased from 0.24, 0.21, and 0.25 m to 0.20, 0.17, and 0.23 m for MP1, MP2, and MP3, respectively. Finally, the average RMSE of the multipath series at the three frequencies was decreased by approximately 19.5%, 20.2%, and 7.5%, respectively. As the multipath is related to the environment around the station, the multipath effects differ among the different stations, which results in the difference between Figure 6 and Table 1.

Multipath Elimination
According to Section 3.2.1, the low-frequency component of the multipath series derived by the wavelet transform can greatly decrease the effect of the multipath error. To evaluate the improvement due to the low-frequency removal, the root mean square error (RMSE) of the one-week multipath series was computed for both without and with the multipath correction, respectively. Figure 6 shows the RMSE of the GEO satellite (C04) at station CUT0. The blue bars represent the RMSE without the multipath error correction, while the red bars represent the RMSE with the multipath error correction. It is obvious that the RMSE of the multipath series considerably decreases with the low-frequency components of the previous day, being removed from the raw time series. The RMSE of the multipath series at station CUT0 on the frequencies of B1, B2, and B3 is decreased by 42.0%, 34.0%, and 60.0% according to the average RMSE, respectively. The RMSE of MP3 is smaller than that of MP1 and MP2, mainly because of its low noise of observable. Table 1 shows the statistical results of the multipath series without and with the multipath correction for each GEO satellite from all of the stations. With the multipath correction, the accuracy of the multipath series at three frequencies is improved. For example, the RMSE of C01 satellite is decreased from 0.24, 0.21, and 0.25 m to 0.20, 0.17, and 0.23 m for MP1, MP2, and MP3, respectively. Finally, the average RMSE of the multipath series at the three frequencies was decreased by approximately 19.5%, 20.2%, and 7.5%, respectively. As the multipath is related to the environment around the station, the multipath effects differ among the different stations, which results in the difference between Figure 6 and Table 1.  Considering that the performance of the multipath correction is related to the time delay from the corrected day, we deeply evaluated the effect of the multipath elimination by using corrections from different time delays. The low-frequency components were extracted from the multipath series with time delays, from one-day to seven-days. Figure 7 shows the statistical results of the C04 satellite in frequencies B1 and B2, and the vertical axis shows the improvement with the multipath correction. From the figure, it is clear that with the time delay increase from one-day to seven-days, the improvement of the multipath RMSE gradually decreases. For example, the best results were obtained using the low-frequency component extracted from the previous day (delay time of one-day). The largest improvement could reach up to 0.3 m, and the average improvements were about 0.05 m and 0.04 m for MP1 and MP2, respectively. However, the results obtained by the corrections from seven-days ago were the worst. The multipath error become even worse for some stations when the corrections were derived from several days ago. This is because with the time extension increase, the correlation of the multipath series between the day correction derived and the day correction applied was decreased. Thus, to eliminate the multipath error as much as possible, we suggest that the latest observable is used to derive the multipath correction.  Considering that the performance of the multipath correction is related to the time delay from the corrected day, we deeply evaluated the effect of the multipath elimination by using corrections from different time delays. The low-frequency components were extracted from the multipath series with time delays, from one-day to seven-days. Figure 7 shows the statistical results of the C04 satellite in frequencies B1 and B2, and the vertical axis shows the improvement with the multipath correction. From the figure, it is clear that with the time delay increase from one-day to seven-days, the improvement of the multipath RMSE gradually decreases. For example, the best results were obtained using the low-frequency component extracted from the previous day (delay time of one-day). The largest improvement could reach up to 0.3 m, and the average improvements were about 0.05 m and 0.04 m for MP1 and MP2, respectively. However, the results obtained by the corrections from seven-days ago were the worst. The multipath error become even worse for some stations when the corrections were derived from several days ago. This is because with the time extension increase, the correlation of the multipath series between the day correction derived and the day correction applied was decreased. Thus, to eliminate the multipath error as much as possible, we suggest that the latest observable is used to derive the multipath correction. Figure 7. Improvement of the multipath series (RMSE) by correcting the multipath errors with different time delays (from one-day to seven-days).

Validation of BDS Precise Data Processing
From Section 3.2, it is validated that the low-frequency component extracted from the multipath series can greatly eliminate the multipath error, and the amount of improvement will decrease with the increase in time span between the day correction derived and the day correction applied. In this section, several experiments were used to validate the performance of the multipath correction, such as the ionospheric delay extraction and satellites clock estimation. All of these experiments were performed in simulated real-time mode, and the time span was chosen as one day. Moreover, they were all performed with two schemes, that is, using the original observables without and with the multipath correction.

Ionospheric Delay Extraction
A geometry-free combination eliminates the geometry-related errors, which are widely used in ionospheric delay estimation [37], and can be described as follows: where , and , denote the noise corresponding to the observables' geometry-free combination of code and phase, respectively. The hardware delay was generally treated as a constant value during a period. When there was no cycle jump in the observables, the combination of ambiguity was also a constant. Ignoring the influence of the high-order terms of the ionospheric delay, the hardware delay and ambiguity can be obtained by , and Φ , , as shown in the following equation.
where 〈•〉 is an average operator for multi-epochs to eliminate noises. The ionospheric delay of the signal propagation path can be obtained by using the observables of the continuous arc through the phase smoothing code range.
where , − , and − are differential hardware delays, which are also named differential code biases (DCBs). The satellite and receiver DCBs can be corrected by the products provided by IGS and averaged through different satellites, respectively.

Validation of BDS Precise Data Processing
From Section 3.2, it is validated that the low-frequency component extracted from the multipath series can greatly eliminate the multipath error, and the amount of improvement will decrease with the increase in time span between the day correction derived and the day correction applied. In this section, several experiments were used to validate the performance of the multipath correction, such as the ionospheric delay extraction and satellites clock estimation. All of these experiments were performed in simulated real-time mode, and the time span was chosen as one day. Moreover, they were all performed with two schemes, that is, using the original observables without and with the multipath correction.

Ionospheric Delay Extraction
A geometry-free combination eliminates the geometry-related errors, which are widely used in ionospheric delay estimation [37], and can be described as follows: where ε P,GF and ε Φ,GF denote the noise corresponding to the observables' geometry-free combination of code and phase, respectively. The hardware delay was generally treated as a constant value during a period. When there was no cycle jump in the observables, the combination of ambiguity was also a constant. Ignoring the influence of the high-order terms of the ionospheric delay, the hardware delay and ambiguity can be obtained by P s r,GF and Φ S r,GF , as shown in the following equation.
where · is an average operator for multi-epochs to eliminate noises. The ionospheric delay of the signal propagation path can be obtained by using the observables of the continuous arc through the phase smoothing code range.
where b r,i − b r,j and b s i − b s j are differential hardware delays, which are also named differential code biases (DCBs). The satellite and receiver DCBs can be corrected by the products provided by IGS and averaged through different satellites, respectively.
According to the above description, the accuracy of the code observables has an impact on the accuracy of the slant path ionospheric delay extraction. The observables of 54 stations within the latitude range of W60 • to E145 • were used in estimating the ionospheric delay in order to study the influences of the multipath error on ionospheric delay extraction, and the observable period is from DOY 46 to 66 in 2017 (16 February 2017 to 7 March 2017). The ionospheric delay estimation is also performed by using the two schemes (i.e., without and with multipath correction). In addition, as the accuracy of the ionospheric delays calculated by PPP mainly depend on phase observables, they are used as the true values in order to evaluate the accuracy of the ionospheric delay obtained by geometry-free combination. Figure 8 shows the time series of the ionospheric residuals on the slant path at station CUT0 without and with the correction of multipath errors. From the figure, it is clear that the accuracy of the ionospheric delay is considerably improved by correcting the multipath error. For example, without the multipath correction, the ionospheric residuals may reach up to eight total electron content unit (TECU) at the initial time of processing. However, with the multipath correction, the ionospheric delay residuals were only approximately three TECU at the initial time. This is because the error of code observables mainly contains white noise, which can be decreased with multi-epochs observables after the low-frequency component of the multipath is removed. Finally, the RMSE of the ionospheric residuals in the first 12-hours of station CUT0 are decreased from 2.54 TECU to 1.46 TECU, which is improved by approximately 42.4%.
According to the above description, the accuracy of the code observables has an impact on the accuracy of the slant path ionospheric delay extraction. The observables of 54 stations within the latitude range of W60° to E145° were used in estimating the ionospheric delay in order to study the influences of the multipath error on ionospheric delay extraction, and the observable period is from DOY 46 to 66 in 2017 (16 February 2017 to 7 March 2017). The ionospheric delay estimation is also performed by using the two schemes (i.e., without and with multipath correction). In addition, as the accuracy of the ionospheric delays calculated by PPP mainly depend on phase observables, they are used as the true values in order to evaluate the accuracy of the ionospheric delay obtained by geometry-free combination. Figure 8 shows the time series of the ionospheric residuals on the slant path at station CUT0 without and with the correction of multipath errors. From the figure, it is clear that the accuracy of the ionospheric delay is considerably improved by correcting the multipath error. For example, without the multipath correction, the ionospheric residuals may reach up to eight total electron content unit (TECU) at the initial time of processing. However, with the multipath correction, the ionospheric delay residuals were only approximately three TECU at the initial time. This is because the error of code observables mainly contains white noise, which can be decreased with multi-epochs observables after the low-frequency component of the multipath is removed. Finally, the RMSE of the ionospheric residuals in the first 12-hours of station CUT0 are decreased from 2.54 TECU to 1.46 TECU, which is improved by approximately 42.4%.

Satellite Clock Estimation
The observables of the 83 stations in the period of DOY 33 to 47 in 2017 (2 February 2017 to 16 February 2017) were used for the BDS satellite clock estimation in order to verify the effect of correcting the GEO's multipath error. The BDS satellite clock was estimated in the simulated real-time mode, and they were set as re-initialization for each day. For the other details of the processing strategy, we refer to the study of [38]. Finally, the multi-system satellite clocks final products provided by Deutsches GeoForschungsZentrum (GFZ) were used to evaluate the accuracy of the satellite clocks' estimation [39]. Figure 9 shows the time series of the BDS satellite clock errors without and with multipath correction. As we were concerned about the convergence time of the satellite clock, the biases of satellite clock between the two products were removed. As shown in Figure 9, with the multipath correction being used, the accuracy of the satellite clock was considerably improved at the first several hours. For example, when the convergence time was approximately two hours, the satellite clock errors of C02 were approximately 0.21 ns and 0.12 ns for without and with multipath correction, respectively. For the satellite clock estimation without multipath correction, 1.5, 2.2, and 3.1 hours were needed in order to converge to 0.2 ns for the C01, C02, and C03 satellites, respectively; however, for the satellite clock estimation with multipath correction, it only needed 0.81, 1.32, and 0.20 hours, respectively. The convergence time was shortened by 46.0%, 40.0%, and 93.5% for the C01, C02, and C03 satellites, respectively. Since C04 and C05 satellites are located in the westernmost and easternmost parts respectively; hence, the convergence time of their satellite clocks was evidently larger than that of the C01, C02, and C03 satellites. Thus, the improvements of the C04 and C05 satellite clock were not evident when the threshold is 0.2 ns. When the threshold was set as 0.5 ns, the convergence times of the C04 and C05 satellite clock were shortened by about 17.6% and 20.2%, respectively.

Satellite Clock Estimation
The observables of the 83 stations in the period of DOY 33 to 47 in 2017 (2 February 2017 to 16 February 2017) were used for the BDS satellite clock estimation in order to verify the effect of correcting the GEO's multipath error. The BDS satellite clock was estimated in the simulated realtime mode, and they were set as re-initialization for each day. For the other details of the processing strategy, we refer to the study of [38]. Finally, the multi-system satellite clocks final products provided by Deutsches GeoForschungsZentrum (GFZ) were used to evaluate the accuracy of the satellite clocks' estimation [39]. Figure 9 shows the time series of the BDS satellite clock errors without and with multipath correction. As we were concerned about the convergence time of the satellite clock, the biases of satellite clock between the two products were removed. As shown in Figure 9, with the multipath correction being used, the accuracy of the satellite clock was considerably improved at the first several hours. For example, when the convergence time was approximately two hours, the satellite clock errors of C02 were approximately 0.21 ns and 0.12 ns for without and with multipath correction, respectively. For the satellite clock estimation without multipath correction, 1.5, 2.2, and 3.1 hours were needed in order to converge to 0.2 ns for the C01, C02, and C03 satellites, respectively; however, for the satellite clock estimation with multipath correction, it only needed 0.81, 1.32, and 0.20 hours, respectively. The convergence time was shortened by 46.0%, 40.0%, and 93.5% for the C01, C02, and C03 satellites, respectively. Since C04 and C05 satellites are located in the westernmost and easternmost parts respectively; hence, the convergence time of their satellite clocks was evidently larger than that of the C01, C02, and C03 satellites. Thus, the improvements of the C04 and C05 satellite clock were not evident when the threshold is 0.2 ns. When the threshold was set as 0.5 ns, the convergence times of the C04 and C05 satellite clock were shortened by about 17.6% and 20.2%, respectively.

Conclusion
The multipath errors cannot be eliminated by the differential technique, because they are closely related to the environment around the station. The existence of multipath errors will greatly decrease the accuracy of the GNSS data processing, especially for the BDS GEO satellites. To correct the multipath error and improve the accuracy of the BDS precise data processing, we deeply investigated

Conclusions
The multipath errors cannot be eliminated by the differential technique, because they are closely related to the environment around the station. The existence of multipath errors will greatly decrease the accuracy of the GNSS data processing, especially for the BDS GEO satellites. To correct the multipath error and improve the accuracy of the BDS precise data processing, we deeply investigated the multipath characteristics of the BDS GEO satellites by using multipath combinations. Given that the multipath combinations consist of noises and low-frequency multipath errors, wavelet transform is used to derive the low-frequency components of the multipath errors. Then, the low-frequency components of the multipath from the first day can be used to correct the multipath errors in real-time.
To validate the method proposed, the multipath series was initially analyzed. The results show that the multipath series was close to white noise when multipath correction was used. The average RMSE of the multipath series was decreased by approximately 19.5%, 20.2%, and 7.5% for B1, B2, and B3, respectively. However, the improvements decreased with the increase in span time delay. For example, the RMSE of the multipath will even become worse for some stations when the time delay is seven-days. Thus, we suggest that the multipath correction should be derived by the latest observables. Moreover, BDS precise data processing, including ionospheric delay extraction and satellite clock estimation, is conducted in simulated real-time mode in order to validate the effectiveness of the multipath correction. For the ionospheric delay extraction, the ionospheric delay calculated by PPP was used as the reference to evaluate the accuracy of the ionospheric delay estimated without and with multipath correction. The results show that the average RMSE of the ionospheric delay on the slant path is decreased by 15.5% with eliminating the multipath error. As for the satellite clock estimation, the results demonstrate that the multipath correction can greatly shorten the convergence time. For example, the times for the satellite clocks were shortened by 46.0%, 40.0%, and 93.5% to converge to 0.2 ns for C01, C02, and C03, respectively. Furthermore, they were shortened by 17.6% and 20.2% to converge to 0.5 ns for the C04 and C05 satellites, respectively. All of these results show that the proposed method can be used to correct multipath error at permanent GNSS stations, which is helpful for real-time BDS precise data processing, such as ionospheric delay estimation and satellite clock estimation. It is also expected that with the multipath error of the IGSO/MEO satellites being investigated, it will improve the accuracy of the BDS precise data processing further.
In order to satisfy the analysis of the complex signals, the theory of time-frequency transformation is further promoted, and wavelet transform is thus produced. The continuous wavelet transform of a signal f is as follows: where ψ(t) is the complex conjugate of ψ(t), and the signal can be reconstructed from the following: where the constant satisfies the following admissibility condition: whereψ (ω) is the Fourier transform of the mother wavelet ψ(t), and ω is the signal frequency.
In practical applications such as signal processing, a finite number of data points are usually given. A discrete version of the wavelet transform is then required, where discrete dilation and translation parameters are used. When considering the computational efficiency, the dyadic a and b values are used, as follows: where j and k are integers. For some particular choices of ψ(t), there exists a corresponding discrete wavelet ψ j,k (t) that has good time-frequency localization properties, such that [40] the following: forms an orthonormal basis for L 2 (R). Using the orthonormal basis, any f (t) ∈ L 2 (R) can be expressed as follows: where the discrete wavelet coefficient α j,k is defined by the following: The wavelet transform defined by Equations (A5)-(A8) is the discrete dyadic wavelet transform. In this discrete wavelet analysis, a signal can be decomposed into its approximations and details. The detail at level j is defined as [21] follows: where Z is the set of positive intergers. The approximation at level J is defined as the sum of the details up to that level, as follows: Then, the signal f(t) can be expressed by the following: From the formula above, we can know that the approximations are related to one another by the following: Equations (A11) and (A12) provide the tree structure of a signal, and also the reconstruction procedure of the original signal. The decomposed approximations and details capture the different frequency bands at different levels, giving information that may not be clearly seen in the original data.