Analysis of the Bias on the Beidou GEO Multipath Combinations

The Beidou navigation satellite system is a very important sensor for positioning in the Asia-Pacific region. The Beidou inclined geosynchronous orbit (IGSO) and medium Earth orbit (MEO) satellites have been analysed in some studies previously conducted by other researchers; this paper seeks to gain more insight regarding the geostationary earth orbit (GEO) satellites. Employing correlation analysis, Fourier transformation and wavelet decomposition, we validate whether there is a systematic bias in their multipath combinations. These biases can be observed clearly in satellites C01, C02 and C04 and have a great correlation with time series instead of elevation, being significantly different from those of the Beidou IGSO and MEO satellites. We propose a correction model to mitigate this bias based on its daily periodicity characteristic. After the model has been applied, the performance of the positioning estimations of the eight stations distributed in the Asia-Pacific region is evaluated and compared. The results show that residuals of multipath series behaves random noise; for the single point positioning (SPP) and precise point positioning (PPP) approaches, the positioning accuracy in the upward direction can be improved by 8 cm and 6 mm, respectively, and by 2 cm and 4 mm, respectively, for the horizontal component.


Introduction
The Beidou navigation system has been providing position, navigation and timing (PNT) services since 27 December 2012, focusing on the whole Asia-Pacific region. The system constellation comprises 15 launched satellites, out of which 13 (five GEO + five IGSO + three MEO) were fully operational in 2015. It has the ability to transmit signals centred on the B2 (1561.098 MHz), B7 (1207.14 MHz) and B6 (1268.52 MHz) signals.
In recent years, several studies and experiments on Beidou have been conducted, particularly on the precise applications [1][2][3][4][5]. Wuhan University can now provide the daily Beidou precise ephemeris, determined using the BeiDou Experimental Tracking Stations [BETS] data, and the orbit precision of the IGSO and MEO satellites can now reach 10 cm in the radial direction [4]. The comparable accuracies prove that as a GPS system Beidou can be a powerful tool for precise positioning. Particularly in the Asia-Pacific region, Beidou is considered as a relevant and valuable complement to legacy systems for establishing improved GNSS services in this area [2]. Current studies indicate that the real-time horizontal component can reach 5 cm, while the vertical component is within 10 cm after the convergence period [6]. Using the PPP approach to determine zenith tropospheric delay (ZTD), the

The Fundamentals of Detecting the Code Bias
The code bias has the appearance of multipath delay, which is always hidden in the multipath series. To detect the code bias, we should obtain the corresponding multipath series first. A so-called multipath combination (MP) is often employed to estimate the multipath delay; the MP delay on frequency i can be expressed as: where P and φ denote the code and phase measurements, respectively; λ is the wavelength; i, j, k are used to denote different frequencies. The linear factor η (i,j,k) is calculated as: η pi,j,kq The linear factor is selected in such a way that the ionospheric delay and non-dispersive delays are cancelled out, leaving the ambiguity differences between the different phase measurements lumped together in B, which is determined as average over raw MP values for each continuous ambiguity block. As a consequence, no absolute MP values are known, just the MP variations over continuous ambiguity blocks. However, long time changes in the pseudorange observations can still be found in the MP time series, which is the basis of our detecting code bias. In most practical applications, setting i = j or i = k can make the computation more convenient.
The theoretical maximum multipath effect on carrier phase is only one quarter of the carrier wave cycle. For Beidou, the maximum multipath effects for B2, B7 and B6 are 4.8, 6.2 and 5.9 cm, respectively. With respect to the magnitude of code bias that is at the level of decimetres or even metres, the multipath effects on the carrier phase are negligible.

The Analysis of the Geo Satellite-Based Code Bias
The multipath bias can be comprehensively analysed using different mathematical tools. Spearman correlation, Fourier transformation and wavelet analyses are all employed in this section. The Spearman correlation coefficient is a distribution-free rank statistic for evaluating the strength of the monotone association between two independent variables. Comparing with other correlation approaches, e.g., Pearson analysis, the Spearman correlation, which operates on the ranks of the data rather than the raw data, does not require a linear relationship between variables. Moreover, Spearman correlation estimation is insensitive to extreme data, which is an advantage because the measurement errors in the pseudorange observables will be neglected by this algorithm [21]. In Section 3.1, we use the Spearman correlation-analysis to determine the correlations between the code bias in Beidou GEO satellites and elevation, as well as code bias and time series. The Fourier transformation is able to depict the multipath time series through amplitude versus frequency, revealing the frequency and the amplitude of each component [22]. In Section 3.2, we take the multipath series as a discrete signal and use the Fourier transformation to investigate the spectral characteristic of multipath time series in the frequency domain, which helps to further identify the periodicity characteristic of the code bias hidden in the multipath series. Wavelet analysis allows us to decompose complicated information contained in a signal into elementary functions associated with different time scales and different frequencies and to reconstruct it with high precision and efficiency [23]. In Section 3.3, wavelet transformation is utilized to isolate the signal with daily periodicity from those at other frequencies; after that, we use the least squares method to estimate the model unknowns. The following sections will explicitly describe the three mentioned approaches.

The Comparison of GEO-Based Code Bias with Other Systems
To better understand the characteristic of the Beidou GEO code bias, we compare the multipath effects induced from GEO satellites with those from GPS and Galileo systems. To make the comparison fair, the Beidou MEO satellites are also selected because all the satellites in the other two navigation systems considered are both MEO satellites. Figure 1 illustrates the multipath series of satellites G04, E11, C12 and C02 and summarizes the correlations between the multipath estimation and the elevation, as well as the time series. From Figure 1a-c, it can be observed that the amplitude and dispersion of the code multipath of the MEO satellites increase significantly as the elevation decreases. The multipath effect nearly increases to 2 m on the first frequency for the GPS system when the elevation is below 20 ¥ and the first frequency observable is more sensitive to the multipath effect than the second frequency. The code multipath of the Galileo satellite has the smallest amplitude and dispersion, for which the multipath variations are almost below 0.5 m and seem independent of elevation. Compared with (a) and (b), the multipath estimation of the Beidou MEO satellite (c) has a trough in all the time series when the satellite has the maximum elevation, and the amplitude becomes bigger when the elevation decreases, which suggests that there is an elevation-dependent bias in the multipath delay. In regard to the Beidou GEO satellite (d), the magnitude of multipath variation is many times bigger than that of the elevation variation, and the trend of multipath series behaves somewhat systematically correlated with time series, which are evidently different from the Beidou MEO satellites. These results agree well with some mentioned references.
Sensors 2016, 16,1252 4 of 22 Figure 1 illustrates the multipath series of satellites G04, E11, C12 and C02 and summarizes the correlations between the multipath estimation and the elevation, as well as the time series. From Figure 1a-c, it can be observed that the amplitude and dispersion of the code multipath of the MEO satellites increase significantly as the elevation decreases. The multipath effect nearly increases to 2 m on the first frequency for the GPS system when the elevation is below 20° and the first frequency observable is more sensitive to the multipath effect than the second frequency. The code multipath of the Galileo satellite has the smallest amplitude and dispersion, for which the multipath variations are almost below 0.5 m and seem independent of elevation. Compared with (a) and (b), the multipath estimation of the Beidou MEO satellite (c) has a trough in all the time series when the satellite has the maximum elevation, and the amplitude becomes bigger when the elevation decreases, which suggests that there is an elevation-dependent bias in the multipath delay. In regard to the Beidou GEO satellite (d), the magnitude of multipath variation is many times bigger than that of the elevation variation, and the trend of multipath series behaves somewhat systematically correlated with time series, which are evidently different from the Beidou MEO satellites. These results agree well with some mentioned references. The code multipath estimations on different frequencies of GPS, Galileo and Beidou satellites show almost the same pattern, while the pattern of the counterparts of GEO satellite from the first frequency has a little deviation with respect to those from the second and the third frequencies. It suggests that the GEO code biases which are absorbed in the multipath delay of the three frequencies need discriminative treatments, taking their different behaviors into account.    The code multipath estimations on different frequencies of GPS, Galileo and Beidou satellites show almost the same pattern, while the pattern of the counterparts of GEO satellite from the first frequency has a little deviation with respect to those from the second and the third frequencies. It suggests that the GEO code biases which are absorbed in the multipath delay of the three frequencies need discriminative treatments, taking their different behaviors into account. Figure 2 shows the Spearman cross-correlations between the multipath delay and the elevation, as well as time series. For Beidou non-MEO satellites, nearly all the correlations with elevation are negatively related, except those for C07 and C12 on the third frequency, and the correlations between multipath delay and time series are quite low and only partially positive, considering different satellites and frequencies. These results, to some extent, verify the existence of elevation-dependent code bias in the multipath delay. For Beidou GEO satellites, the correlation results show that the multipath delay has correlations with both time series and elevations. However, the GEO satellites show daily variations in their fixed frames of the order of approximately 2000 km, which results in their elevation also having the daily variation periodicity characteristic [11]. The variation of GEO code bias should be carefully analysed to be reasonably mitigated. For this purpose, it becomes necessary to analyse observations from some stations with different elevations. Although some conclusions have been made before [11,17], to guarantee the correctness and reliability of the following findings, we conduct this work again by using newer and more data collected from stations that are distributed over the entire Asia-Pacific region. Based on this, we select 8 stations distributed over the Asia-pacific region as experiment subjects; their distributions are depicted in Figure 3. The observation sampling interval is 30 s. Most of the stations come from the MGEX network, except station BJF1, which comes from the Chinese iGMAS network. To explicitly describe the periodical variation trends of the GEO-based code bias over consecutive days, we select the multipath series from 50 to 55 days of year (Doy), 2015, as representative results to illustrate the multipath effects.
Sensors 2016, 16,1252 5 of 22 Figure 2 shows the Spearman cross-correlations between the multipath delay and the elevation, as well as time series. For Beidou non-MEO satellites, nearly all the correlations with elevation are negatively related, except those for C07 and C12 on the third frequency, and the correlations between multipath delay and time series are quite low and only partially positive, considering different satellites and frequencies. These results, to some extent, verify the existence of elevation-dependent code bias in the multipath delay. For Beidou GEO satellites, the correlation results show that the multipath delay has correlations with both time series and elevations. However, the GEO satellites show daily variations in their fixed frames of the order of approximately 2000 km, which results in their elevation also having the daily variation periodicity characteristic [11]. The variation of GEO code bias should be carefully analysed to be reasonably mitigated. For this purpose, it becomes necessary to analyse observations from some stations with different elevations. Although some conclusions have been made before [11,17], to guarantee the correctness and reliability of the following findings, we conduct this work again by using newer and more data collected from stations that are distributed over the entire Asia-Pacific region. Based on this, we select 8 stations distributed over the Asia-pacific region as experiment subjects; their distributions are depicted in Figure 3. The observation sampling interval is 30 s. Most of the stations come from the MGEX network, except station BJF1, which comes from the Chinese iGMAS network. To explicitly describe the periodical variation trends of the GEO-based code bias over consecutive days, we select the multipath series from 50 to 55 days of year (Doy), 2015, as representative results to illustrate the multipath effects.
The code observations of the first frequency have the best signal quality, being able to clearly and directly present the required results [1,19]. On the contrary, many unrelated high-frequency signals are included in the third frequency observations; the following result (i.e., Table 1) also suggests that the amplitude of periodical variations of the third frequency is the smallest. The performance of the code observations of the second frequency is in between those of the other two frequencies. To show the periodicity characteristic of the code bias, as well as present how to detect and remove unrelated signals, we set the observations of the second frequency as an example to conduct the relevant analysis in Section 3. For simplicity sake, the multipath results of the other frequencies are not displayed here.    Many primary investigations have verified that the satellite-induced code bias is independent of receiver types and free from observation environments [14,18,20]. Moreover, Zhao et al. found that the GEO phase bias, which varied with satellite elevation at the Northern Hemisphere stations, had a trend opposite to that in the Southern Hemisphere [17]. Based on this finding, the pseudorange bias may also have a similar characteristic. Thus, in the following analyses, we select stations JFNG and CUT0, which, respectively represent stations in the Northern and the Southern Hemisphere, to investigate the effects of azimuth and elevation on the code bias.
The multipath estimations of the two stations are presented in Figure 4; the figure shows two phenomena: (1) the elevation variation of the GEO satellites at station JFNG is opposite to that at station CUT0; (2) the variation trends of the GEO multipath series are nearly the same at the two stations. The GEO satellites move around the Equator. When they move closer to one station, the elevation will decrease for the station in the opposite hemisphere, which is the reason for the first phenomenon. The second phenomenon shows that the code bias is independent of the azimuth and the elevation and that the systematic multipath delay is most related to the satellite due to the two stations having obviously different observation environments. Based on these two phenomena, the code bias varying with satellite elevation in the Northern Hemisphere has a trend opposite to that in the Southern Hemisphere, which agrees well with the experimental results regarding the phase bias [19]. The code observations of the first frequency have the best signal quality, being able to clearly and directly present the required results [1,19]. On the contrary, many unrelated high-frequency signals are included in the third frequency observations; the following result (i.e., Table 1) also suggests that the amplitude of periodical variations of the third frequency is the smallest. The performance of the code observations of the second frequency is in between those of the other two frequencies. To show the periodicity characteristic of the code bias, as well as present how to detect and remove unrelated signals, we set the observations of the second frequency as an example to conduct the relevant analysis in Section 3. For simplicity sake, the multipath results of the other frequencies are not displayed here. Many primary investigations have verified that the satellite-induced code bias is independent of receiver types and free from observation environments [14,18,20]. Moreover, Zhao et al. found that the GEO phase bias, which varied with satellite elevation at the Northern Hemisphere stations, had a trend opposite to that in the Southern Hemisphere [17]. Based on this finding, the pseudorange bias may also have a similar characteristic. Thus, in the following analyses, we select stations JFNG and CUT0, which, respectively represent stations in the Northern and the Southern Hemisphere, to investigate the effects of azimuth and elevation on the code bias.
The multipath estimations of the two stations are presented in Figure 4; the figure shows two phenomena: (1) the elevation variation of the GEO satellites at station JFNG is opposite to that at (2) the variation trends of the GEO multipath series are nearly the same at the two stations. The GEO satellites move around the Equator. When they move closer to one station, the elevation will decrease for the station in the opposite hemisphere, which is the reason for the first phenomenon. The second phenomenon shows that the code bias is independent of the azimuth and the elevation and that the systematic multipath delay is most related to the satellite due to the two stations having obviously different observation environments. Based on these two phenomena, the code bias varying with satellite elevation in the Northern Hemisphere has a trend opposite to that in the Southern Hemisphere, which agrees well with the experimental results regarding the phase bias [19].  From the results in Figures 4 and 5, evident daily periodicity can be observed in the multipath time series of C01, C02 and C04. Their initial values and variation trends are similar across the three satellites, further indicating that the time-dependent bias is most related to the satellite instead of the elevation or station. Take stations CUTO and JFNG as an example; these two stations are distantly separated and their observation environments are independent of each other; their elevation trends are even opposite, while the trends of the multipath estimations at the two stations are much similar for each satellite of the three. In this sense, the multipath delays from C01, C02 and C04 might not only contain effects from the surroundings of the station but also contain some bias from the satellite, e.g., internal hardware delay [19]. In Figure 4, the trend of C04 at station CUTO has an obvious characteristic of periodicity variation, but this is not the same for station JFNG; the reason for this phenomenon will be discussed in the following sections. In Figure 5, for the three GEO satellites (C01,C02 and C04), the shape of each multipath series at the 6 stations seems sinusoidal, the shape of each multipath series at the experimental stations seems sinusoidal, which has also been demonstrated in reference [11], and their maximum variation ranges are all limited to 0.12 m to 0.16 m. In this sense, the multipath biases at the 6 stations seem to have similar variation shapes. Figure 6 shows that the variations of C03 and C05 illustrate the random distribution characteristic, which is evidently different from that in the other three GEO satellites.  From the results in Figures 4 and 5, evident daily periodicity can be observed in the multipath time series of C01, C02 and C04. Their initial values and variation trends are similar across the three satellites, further indicating that the time-dependent bias is most related to the satellite instead of the elevation or station. Take stations CUTO and JFNG as an example; these two stations are distantly separated and their observation environments are independent of each other; their elevation trends are even opposite, while the trends of the multipath estimations at the two stations are much similar for each satellite of the three. In this sense, the multipath delays from C01, C02 and C04 might not only contain effects from the surroundings of the station but also contain some bias from the satellite, e.g., internal hardware delay [19]. In Figure 4, the trend of C04 at station CUTO has an obvious characteristic of periodicity variation, but this is not the same for station JFNG; the reason for this phenomenon will be discussed in the following sections. In Figure 5, for the three GEO satellites (C01,C02 and C04), the shape of each multipath series at the 6 stations seems sinusoidal, the shape of each multipath series at the experimental stations seems sinusoidal, which has also been demonstrated in reference [11], and their maximum variation ranges are all limited to 0.12 m to 0.16 m. In this sense, the multipath biases at the 6 stations seem to have similar variation shapes. Figure 6 shows that the variations of C03 and C05 illustrate the random distribution characteristic, which is evidently different from that in the other three GEO satellites.    Figure 6. The multipath series of C03 and C05 at the other six stations. Green is the multipath estimation; red is the Fourier fitting results of the multipath series; blue is the wavelet reconstructed multipath series; the cycle is the edge of reconstructed multipath series between different stations. Figure 7 presents the Spearman cross-correlations between the multipath estimations at station CUT0 and those at station JFNG. The magnitudes of the correlations of C01, C02 and C04 are many times bigger than those of C03 and C05. This phenomenon verifies that there are some biases in satellitesC01 and C02, which affect the pseudorange observations at the two stations in terms of communal errors. C03 and C05 do not show such systematic bias. Because the data of satellite hardware is not available, the possible reason may be that satellites C03 and C05 do have not such bias or these biases are constants, which are absorbed in B in Equation (1). Reference [1] reported that satellite C03 had the best code data quality, with the smallest standard deviation compared with that of C01, C02 and C04 (C05 was not available in that paper); the reason may just be that C03 has no such systematic code bias. The overall magnitudes of all the cross-correlations are not as large as expected, being mainly affected by the code observation noise, whose standard deviations are even larger than the variations of the code bias under some conditions. Figure 7 presents the Spearman cross-correlations between the multipath estimations at station CUT0 and those at station JFNG. The magnitudes of the correlations of C01, C02 and C04 are many times bigger than those of C03 and C05. This phenomenon verifies that there are some biases in satellitesC01 and C02, which affect the pseudorange observations at the two stations in terms of communal errors. C03 and C05 do not show such systematic bias. Because the data of satellite hardware is not available, the possible reason may be that satellites C03 and C05 do have not such bias or these biases are constants, which are absorbed in B in Equation (1). Reference [1] reported that satellite C03 had the best code data quality, with the smallest standard deviation compared with that of C01, C02 and C04 (C05 was not available in that paper); the reason may just be that C03 has no such systematic code bias. The overall magnitudes of all the cross-correlations are not as large as expected, being mainly affected by the code observation noise, whose standard deviations are even larger than the variations of the code bias under some conditions. To quantitatively verify that the biases at the 6 stations are similar for certain satellites, Figures 8 and 9 demonstrate the Spearman cross-correlations between multipath bias at the six stations for all five GEO satellites. Figure 8 presents the multipath correlations for satellites C01, C02 and C04. Here, all the Spearman correlations between multipath series at different stations are positively related. Most of their correlation values are larger than 0.5, except the results of station MRO1, which range from 0.3 to 0.5. These results demonstrate that their multipath biases are highly correlated, which means that the multipath variation trends at the 6 stations will be quite similar; otherwise their correlations will not be positive with such large values. On the contrary, Figure 9 shows that the multipath correlations of C03 or C05 are only partially positive; most correlation values are less than 0.3, which means that the bias variations at the six stations are not similar for satellites C03 and C05. These deductions from Figures 8 and 9 further verify that bias with the daily periodicity characteristic is induced for satellites C01, C02 and C04, while the bias of satellites C03 and C05 has no such characteristic. To quantitatively verify that the biases at the 6 stations are similar for certain satellites, Figures 8 and 9 demonstrate the Spearman cross-correlations between multipath bias at the six stations for all five GEO satellites. Figure 8 presents the multipath correlations for satellites C01, C02 and C04. Here, all the Spearman correlations between multipath series at different stations are positively related. Most of their correlation values are larger than 0.5, except the results of station MRO1, which range from 0.3 to 0.5. These results demonstrate that their multipath biases are highly correlated, which means that the multipath variation trends at the 6 stations will be quite similar; otherwise their correlations will not be positive with such large values. On the contrary, Figure 9 shows that the multipath correlations of C03 or C05 are only partially positive; most correlation values are less than 0.3, which means that the bias variations at the six stations are not similar for satellites C03 and C05. These deductions from Figures 8 and 9 further verify that bias with the daily periodicity characteristic is induced for satellites C01, C02 and C04, while the bias of satellites C03 and C05 has no such characteristic.  Figure 7 presents the Spearman cross-correlations between the multipath estimations at station CUT0 and those at station JFNG. The magnitudes of the correlations of C01, C02 and C04 are many times bigger than those of C03 and C05. This phenomenon verifies that there are some biases in satellitesC01 and C02, which affect the pseudorange observations at the two stations in terms of communal errors. C03 and C05 do not show such systematic bias. Because the data of satellite hardware is not available, the possible reason may be that satellites C03 and C05 do have not such bias or these biases are constants, which are absorbed in B in Equation (1). Reference [1] reported that satellite C03 had the best code data quality, with the smallest standard deviation compared with that of C01, C02 and C04 (C05 was not available in that paper); the reason may just be that C03 has no such systematic code bias. The overall magnitudes of all the cross-correlations are not as large as expected, being mainly affected by the code observation noise, whose standard deviations are even larger than the variations of the code bias under some conditions. To quantitatively verify that the biases at the 6 stations are similar for certain satellites, Figures 8 and 9 demonstrate the Spearman cross-correlations between multipath bias at the six stations for all five GEO satellites. Figure 8 presents the multipath correlations for satellites C01, C02 and C04. Here, all the Spearman correlations between multipath series at different stations are positively related. Most of their correlation values are larger than 0.5, except the results of station MRO1, which range from 0.3 to 0.5. These results demonstrate that their multipath biases are highly correlated, which means that the multipath variation trends at the 6 stations will be quite similar; otherwise their correlations will not be positive with such large values. On the contrary, Figure 9 shows that the multipath correlations of C03 or C05 are only partially positive; most correlation values are less than 0.3, which means that the bias variations at the six stations are not similar for satellites C03 and C05. These deductions from Figures 8 and 9 further verify that bias with the daily periodicity characteristic is induced for satellites C01, C02 and C04, while the bias of satellites C03 and C05 has no such characteristic.

The Spectrum Analysis of GEO Code Bias
Figures 10 and 11 illustrate the multipath results of stations CUT0 and JFNG, as well as their corresponding fitting curves and magnitude spectrums derived from Fourier transformation. Since we will take the observations at station JFNG as example to illustrate how to use waveletreconstruction to reconstruct the raw signal, we resample the data at station JFNG with the interval of 240 s for convenient computation. The sampling interval at the other stations is 30 s. In Figure 10, the fitting curve and magnitude spectrum both illustrate the code bias of C01, C02 and C04 have evident daily periodicity characteristic and similar variation trends. The same holds for station JFNG

The Spectrum Analysis of GEO Code Bias
Figures 10 and 11 illustrate the multipath results of stations CUT0 and JFNG, as well as their corresponding fitting curves and magnitude spectrums derived from Fourier transformation. Since we will take the observations at station JFNG as example to illustrate how to use waveletreconstruction to reconstruct the raw signal, we resample the data at station JFNG with the interval of 240 s for convenient computation. The sampling interval at the other stations is 30 s. In Figure 10, the fitting curve and magnitude spectrum both illustrate the code bias of C01, C02 and C04 have evident daily periodicity characteristic and similar variation trends. The same holds for station JFNG in Figure 4, except satellite C04. The bottom panel in Figure 11 illustrates there is an additional  Figures 10 and 11 illustrate the multipath results of stations CUT0 and JFNG, as well as their corresponding fitting curves and magnitude spectrums derived from Fourier transformation. Since we will take the observations at station JFNG as example to illustrate how to use wavelet-reconstruction to reconstruct the raw signal, we resample the data at station JFNG with the interval of 240 s for convenient computation. The sampling interval at the other stations is 30 s. In Figure 10, the fitting curve and magnitude spectrum both illustrate the code bias of C01, C02 and C04 have evident daily periodicity characteristic and similar variation trends. The same holds for station JFNG in Figure 4, except satellite C04. The bottom panel in Figure 11 illustrates there is an additional semi-day signal also included in the multipath delays of C04. However, during the same experiment time span, other stations e.g., station GMSD ( Figure 12) illustrates the characteristic of daily periodicity, which reveals an excellent agreement with the results about CUT0. It can be deducted that the semi-day signal of C04 come or derived from its own observation environment, and it is most related to satellite. time span, other stations e.g., station GMSD ( Figure 12) illustrates the characteristic of daily periodicity, which reveals an excellent agreement with the results about CUT0. It can be deducted that the semi-day signal of C04 come or derived from its own observation environment, and it is most related to satellite. time span, other stations e.g., station GMSD ( Figure 12) illustrates the characteristic of daily periodicity, which reveals an excellent agreement with the results about CUT0. It can be deducted that the semi-day signal of C04 come or derived from its own observation environment, and it is most related to satellite. In practical application, a series can be reliably believed to have the periodicity characteristic when its corresponding amplitude is more than 80. However, Figure 13 shows that the multipath estimations of C03 and C05 could not satisfy this condition at the two stations at the same time, to a great extent, which further suggests that there is no periodicity bias in satellites C03 and C05. In practical application, a series can be reliably believed to have the periodicity characteristic when its corresponding amplitude is more than 80. However, Figure 13 shows that the multipath estimations of C03 and C05 could not satisfy this condition at the two stations at the same time, to a great extent, which further suggests that there is no periodicity bias in satellites C03 and C05. In practical application, a series can be reliably believed to have the periodicity characteristic when its corresponding amplitude is more than 80. However, Figure 13 shows that the multipath estimations of C03 and C05 could not satisfy this condition at the two stations at the same time, to a great extent, which further suggests that there is no periodicity bias in satellites C03 and C05. In practical application, a series can be reliably believed to have the periodicity characteristic when its corresponding amplitude is more than 80. However, Figure 13 shows that the multipath estimations of C03 and C05 could not satisfy this condition at the two stations at the same time, to a great extent, which further suggests that there is no periodicity bias in satellites C03 and C05. The above results show that C01, C02 and C04 mostly have satellite-induced code bias with a clear periodicity variation characteristic. Satellites C03 and C05 have no such bias; perhaps the magnitude of such bias is too small and can be ignored or their systematic bias is constant, being absorbed by term B in Equation (1).

The Spectrum Analysis of GEO Code Bias
Assume the sampling interval as Δ (ins), frequency as f and periodicity as T. For the multipath series at stations CUT0 and GMSD, the Δ is 30 s, while f, which corresponds to the first peak (also the most evident peak), is approximately 3.47 × 10 -4 (Hz). For the multipath series at stations JFNG, the Δ is 240 s, while f is approximately 2.78 × 10 -3 (Hz). Using Equation (3), all the periodicities of the multipath combinations of satellites C01, C02 and C04 are nearly equal to 86,400 s, which further quantitatively verifies that the periodicity of their corresponding multipath bias is one day: Moreover, a semi-day signal could be detected in all the mentioned magnitude spectrums. Although their corresponding amplitudes are much smaller than those of the signals with the daily periodicity characteristic, such semi-day signal is nearly observable in all the experimental stations. Because the data of the satellite internal hardware is not available, we can only assume that the semiday signal is most related to satellites.

Time Series-Depending Model for GEO Code Bias
Based on the above analysis, the daily repeating code bias plays an important role in satellites C01, C02 and C04. It is necessary and significant to model this time-dependent bias to mitigate its effect on the positioning. The Beidou observations collected from 040 to 060 (Doy), 2015, at eight stations are used; we resample the data at the interval of 240 s for convenient computation. The computation process can be concluded as follows: (1) Classify all the observables according to the types of frequency and satellite; after that, we use the wavelet transformation to remove the signal trends; (2) There have many unrelated high-frequency signals in the raw observations, which disturb the analysis of the special signals that have the daily periodicity characteristic. Considering this, we use wavelet-reconstruction to reconstruct the raw signals. Because the sampling interval (Δ) of the experimental data is 240 s, we use the Daubechies wavelet with a 45 degree vanishing moment to analyse and construct the multipath series. Thus, we obtain the signal with a frequency band of 8 9 2~2  , i.e., 17.06~34.12 h.
After step 2, the reconstructed signal can clearly illustrate the characteristic of daily periodicity. For example, the multipath series of C04 at station JFNG does not illustrate an obvious daily periodicity characteristic in Figure 4 due to a semi-daily signal being included, which is also mentioned in Figure 11. After the second step, Figure 14 shows that the daily periodicity characteristic can clearly be observed in the reconstructed signals. The above results show that C01, C02 and C04 mostly have satellite-induced code bias with a clear periodicity variation characteristic. Satellites C03 and C05 have no such bias; perhaps the magnitude of such bias is too small and can be ignored or their systematic bias is constant, being absorbed by term B in Equation (1).
Assume the sampling interval as ∆ (ins), frequency as f and periodicity as T. For the multipath series at stations CUT0 and GMSD, the ∆ is 30 s, while f, which corresponds to the first peak (also the most evident peak), is approximately 3.47 ¢ 10 ¡4 (Hz). For the multipath series at stations JFNG, the ∆ is 240 s, while f is approximately 2.78 ¢ 10 ¡3 (Hz). Using Equation (3), all the periodicities of the multipath combinations of satellites C01, C02 and C04 are nearly equal to 86,400 s, which further quantitatively verifies that the periodicity of their corresponding multipath bias is one day: Moreover, a semi-day signal could be detected in all the mentioned magnitude spectrums. Although their corresponding amplitudes are much smaller than those of the signals with the daily periodicity characteristic, such semi-day signal is nearly observable in all the experimental stations. Because the data of the satellite internal hardware is not available, we can only assume that the semi-day signal is most related to satellites.

Time Series-Depending Model for GEO Code Bias
Based on the above analysis, the daily repeating code bias plays an important role in satellites C01, C02 and C04. It is necessary and significant to model this time-dependent bias to mitigate its effect on the positioning. The Beidou observations collected from 040 to 060 (Doy), 2015, at eight stations are used; we resample the data at the interval of 240 s for convenient computation. The computation process can be concluded as follows: (1) Classify all the observables according to the types of frequency and satellite; after that, we use the wavelet transformation to remove the signal trends; (2) There have many unrelated high-frequency signals in the raw observations, which disturb the analysis of the special signals that have the daily periodicity characteristic. Considering this, we use wavelet-reconstruction to reconstruct the raw signals. Because the sampling interval (∆) of the experimental data is 240 s, we use the Daubechies wavelet with a 45 degree vanishing moment to analyse and construct the multipath series. Thus, we obtain the signal with a frequency band of 2 8 ∆ 2 9 ∆, i.e., 17.06~34.12 h.
After step 2, the reconstructed signal can clearly illustrate the characteristic of daily periodicity. For example, the multipath series of C04 at station JFNG does not illustrate an obvious daily periodicity characteristic in Figure 4 due to a semi-daily signal being included, which is also mentioned in Figure 11.
After the second step, Figure 14 shows that the daily periodicity characteristic can clearly be observed in the reconstructed signals.
(3) After many fitting methods were tested, we found that the trigonometric function is the most effective for this study because this method yields the smallest RMS value among the tested fitting approaches. The estimation of code bias Cb can be modelled as: where x represents the time difference between observation time and midnight on 050 (Doy), 2015, with the unit of seconds, and ω, a, b are the unknown parameters. To facilitate the least squares adjustment, we transform Equation (4) into the standard indirect adjustment model; it is expressed as:  Table 1. The results indicate that the magnitudes of bias corrections of the first and second frequencies are at the same level, both being bigger than that of the third frequency. The periodicity signal, which is extracted from the multipath series of the third frequency, has the lowest amplitude. We use the equation: to estimate periodicity T. The results show that periodicity varies from 86,069 to 86,237 s for different frequencies and satellites, being roughly equal to the length of one sidereal day, i.e., 86,160 s. Note that the estimations of ω, a, b will change with the time varies; details of the relevant deductions will be discussed in another paper. (3) After many fitting methods were tested, we found that the trigonometric function is the most effective for this study because this method yields the smallest RMS value among the tested fitting approaches. The estimation of code bias Cb can be modelled as:

Experiments
Turning now to the evaluation of the performance of the GEO code bias model, single point positioning (SPP) and precise point positioning (PPP), with 20 days of Beidou observations, are conducted and the results are compared. We mainly describe the results at station CUT0 for simplicity. During the experiments, the precise station coordinates come from weekly solutions, and positioning errors can be computed by subtracting the estimated coordinates from the precisely estimated coordinates [24,25].
For comparison with the multipath series of satellites C01, C02 and C04 in Figure 4, the right panel in Figure 15 depicts the multipath estimations to which the time-dependent correction model has been applied. After applying the model, the time-dependent periodicity characteristic in the previous multipath estimations almost disappears. The standard deviation (std) decreases from 0.29 m to 0.22 m for satellite C01 and decreases by 0.06 m and 0.09 m for C02 and C04, respectively. It can be deduced that the signal quality is improved after the model has been applied. For comparison with the multipath series of satellites C01, C02 and C04 in Figure 4, the right panel in Figure 15 depicts the multipath estimations to which the time-dependent correction model has been applied. After applying the model, the time-dependent periodicity characteristic in the previous multipath estimations almost disappears. The standard deviation (std) decreases from 0.29 m to 0.22 m for satellite C01 and decreases by 0.06 m and 0.09 m for C02 and C04, respectively. It can be deduced that the signal quality is improved after the model has been applied.
SPP, which uses code measurements to produce coordinates, is usually affected by the code bias. Figures 16 and 17 show that the correction model has certain effects on the SPP estimations. After the model is applied, the positioning bias decreases by 0.13 m in the U direction and by approximately 0.03 m in the horizontal components. A decrease of more than 0.05 m is also found in terms of the standard deviation in the upward (U) direction. The correction model performs better regarding the improvement of position precision in the U direction rather than in the horizontal components. SPP, which uses code measurements to produce coordinates, is usually affected by the code bias. Figures 16 and 17 show that the correction model has certain effects on the SPP estimations. After the model is applied, the positioning bias decreases by 0.13 m in the U direction and by approximately 0.03 m in the horizontal components. A decrease of more than 0.05 m is also found in terms of the standard deviation in the upward (U) direction. The correction model performs better regarding the improvement of position precision in the U direction rather than in the horizontal components.   Figure 18 summarizes the average results of the 20 days for each experimental station; it is based on ionosphere-free combination of the first and second frequencies. The improvement of positioning precision in the U direction is the highest, and the precision improvement in the east (E) direction is lower than that in the U direction. The precision improvement in the north (N) direction is the lowest, which may be caused by the geometric structure of the three corrected GEO satellites. The overall correction values in Figure 19 are smaller than those in Figure 18, possibly because the correction value of the third frequency is smaller than that of the second frequency.    Figure 18 summarizes the average results of the 20 days for each experimental station; it is based on ionosphere-free combination of the first and second frequencies. The improvement of positioning precision in the U direction is the highest, and the precision improvement in the east (E) direction is lower than that in the U direction. The precision improvement in the north (N) direction is the lowest, which may be caused by the geometric structure of the three corrected GEO satellites. The overall correction values in Figure 19 are smaller than those in Figure 18, possibly because the correction value of the third frequency is smaller than that of the second frequency.  Figure 17. "Positioning difference" is the difference between SPP estimations with and without employing the correction model. The "mean bias" represents the equal result of the positioning errors. "Uncorrected" indicates the results when the code bias is not corrected; "Corrected" indicates the results after the model has been applied. Figure 18 summarizes the average results of the 20 days for each experimental station; it is based on ionosphere-free combination of the first and second frequencies. The improvement of positioning precision in the U direction is the highest, and the precision improvement in the east (E) direction is lower than that in the U direction. The precision improvement in the north (N) direction is the lowest, which may be caused by the geometric structure of the three corrected GEO satellites. The overall correction values in Figure 19 are smaller than those in Figure 18, possibly because the correction value of the third frequency is smaller than that of the second frequency. The pseudorange measurements occupy 1/10000 weights in the static PPP approach with respect to the phase observables, which makes the pseudorange observations contribute marginally to the final positioning precision [7]. However, the pseudorange has a great effect on the convergence efficiency; Figure 20 illustrates the positioning errors during the first 180 epochs at station CUTO. After converging with 80 epochs, the positioning precision, which is based on corrected pseudorange observations, improves by 0.16 m, 0.08 m and 0.01 m with respect to that using the raw observables in the U, E and N directions, respectively. The pseudorange measurements occupy 1/10000 weights in the static PPP approach with respect to the phase observables, which makes the pseudorange observations contribute marginally to the final positioning precision [7]. However, the pseudorange has a great effect on the convergence efficiency; Figure 20 illustrates the positioning errors during the first 180 epochs at station CUTO. After converging with 80 epochs, the positioning precision, which is based on corrected pseudorange observations, improves by 0.16 m, 0.08 m and 0.01 m with respect to that using the raw observables in the U, E and N directions, respectively. The pseudorange measurements occupy 1/10000 weights in the static PPP approach with respect to the phase observables, which makes the pseudorange observations contribute marginally to the final positioning precision [7]. However, the pseudorange has a great effect on the convergence efficiency; Figure 20 illustrates the positioning errors during the first 180 epochs at station CUTO. After converging with 80 epochs, the positioning precision, which is based on corrected pseudorange observations, improves by 0.16 m, 0.08 m and 0.01 m with respect to that using the raw observables in the U, E and N directions, respectively.     The model has a better performance for the combination of the first and second frequencies than for the combination of the first and third frequencies. The improvements of positioning precision in the U direction can be easily observed in the two figures, and the horizontal components have a little deviation and do not seem to have any cross-correlations between the two ionosphere-free combinations. The equal improvements are approximately 6 mm in the U direction and 3 mm for the horizontal components. However, Figure 24 shows one exception, i.e., the model decreases the positioning precision in the N direction at station MAJU. One of the possible explanations can be corrections of the first frequency are not totally in agreement with the truth; the corrections of the second frequency could make up for this shortcoming, while the corrections of the third frequency cannot do so due to values with smaller magnitude. This suggests that these exceptions mainly occur for the combinations of the first and third frequencies and not for the combinations of the first and second frequencies. The model has a better performance for the combination of the first and second frequencies than for the combination of the first and third frequencies. The improvements of positioning precision in the U direction can be easily observed in the two figures, and the horizontal components have a little deviation and do not seem to have any cross-correlations between the two ionosphere-free combinations. The equal improvements are approximately 6 mm in the U direction and 3 mm for the horizontal components. However, Figure 24 shows one exception, i.e., the model decreases the positioning precision in the N direction at station MAJU. One of the possible explanations can be corrections of the first frequency are not totally in agreement with the truth; the corrections of the second frequency could make up for this shortcoming, while the corrections of the third frequency cannot do so due to values with smaller magnitude. This suggests that these exceptions mainly occur for the combinations of the first and third frequencies and not for the combinations of the first and second frequencies.
Sensors 2016, 16,1252 19 of 22 Figure 22. "Positioning difference" is the differences between kinematic PPP estimations with and without employing the correction model. "Uncorrected" indicates the results when the code bias is not corrected; "Corrected" indicates the results after the model has been applied. Figures 23 and 24 illustrate the performance of the model in the kinematic PPP approach. The model has a better performance for the combination of the first and second frequencies than for the combination of the first and third frequencies. The improvements of positioning precision in the U direction can be easily observed in the two figures, and the horizontal components have a little deviation and do not seem to have any cross-correlations between the two ionosphere-free combinations. The equal improvements are approximately 6 mm in the U direction and 3 mm for the horizontal components. However, Figure 24 shows one exception, i.e., the model decreases the positioning precision in the N direction at station MAJU. One of the possible explanations can be corrections of the first frequency are not totally in agreement with the truth; the corrections of the second frequency could make up for this shortcoming, while the corrections of the third frequency cannot do so due to values with smaller magnitude. This suggests that these exceptions mainly occur for the combinations of the first and third frequencies and not for the combinations of the first and second frequencies.

Conclusions
In this paper, the most recent GNSS observables collected from eight stations in the Asia-Pacific region are employed to investigate the characteristic of code bias of the Beidou GEO satellites. Employing some algorithms, such as Spearman correlation analysis, Fourier transformation and wavelet decomposition, this paper validates that the code multipath combination of the BDS satellites not only contains multipath effects from the surroundings of the stations but also some multipathlike biases from the satellites. The multipath bias shows the daily periodical characteristic as well as similar variation trends of the eight experimental stations for GEO satellites C01, C02 and C04; this was not detected in C03 and C05. Experiments also verify that the periodical variation of the multipath bias on satellites C01, C02 and C04 is most related to the satellites rather than the external observation environments. This finding could help verify the hardware questions for GEO satellites C01/C02/C04 through comparison with C03 and C05. By analysing the spectrum amplitude results, we also detected a semi-day signal in our experiments; this signal is observable in nearly all the multipath series. We think that it is an issue that deserves further analysis.
We present the correction model by extracting the multipath combinations from many previous days. Validation experiments verify, after the model has been applied, that the multipath residuals become random and the periodicity characteristic in the previous multipath series is nearly removed. The convergence efficiency is thus improved in the static PPP experiment. Moreover, the precisions of single point positioning and kinematic precise positioning estimations are also improved, especially in the U direction.
In general, this manuscript particularly researched the Beidou GEO satellites, aiming to reveal more characteristics of the code bias in GEO satellites. Based on our findings, more mathematical models and physical settings are expected to be presented to mitigate the GEO multipath delay, thereby improving the Beidou positioning accuracy.

Conclusions
In this paper, the most recent GNSS observables collected from eight stations in the Asia-Pacific region are employed to investigate the characteristic of code bias of the Beidou GEO satellites. Employing some algorithms, such as Spearman correlation analysis, Fourier transformation and wavelet decomposition, this paper validates that the code multipath combination of the BDS satellites not only contains multipath effects from the surroundings of the stations but also some multipath-like biases from the satellites. The multipath bias shows the daily periodical characteristic as well as similar variation trends of the eight experimental stations for GEO satellites C01, C02 and C04; this was not detected in C03 and C05. Experiments also verify that the periodical variation of the multipath bias on satellites C01, C02 and C04 is most related to the satellites rather than the external observation environments. This finding could help verify the hardware questions for GEO satellites C01/C02/C04 through comparison with C03 and C05. By analysing the spectrum amplitude results, we also detected a semi-day signal in our experiments; this signal is observable in nearly all the multipath series. We think that it is an issue that deserves further analysis.
We present the correction model by extracting the multipath combinations from many previous days. Validation experiments verify, after the model has been applied, that the multipath residuals become random and the periodicity characteristic in the previous multipath series is nearly removed. The convergence efficiency is thus improved in the static PPP experiment. Moreover, the precisions of single point positioning and kinematic precise positioning estimations are also improved, especially in the U direction.
In general, this manuscript particularly researched the Beidou GEO satellites, aiming to reveal more characteristics of the code bias in GEO satellites. Based on our findings, more mathematical models and physical settings are expected to be presented to mitigate the GEO multipath delay, thereby improving the Beidou positioning accuracy.