sensors Article

Mitigating BeiDou Satellite-Induced Code Bias: Taking into Account the Stochastic Model of Corrections Fei Guo 1,2,3 , Xin Li 1 and Wanke Liu 1,2,3, * 1 2 3

*

School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China; [email protected] (F.G.); [email protected] (X.L.) Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan 430079, China Collaborative Innovation Center for Geospatial Technology, Wuhan 430079, China Correspondence: [email protected]; Tel.: +86-27-6875-8530

Academic Editor: Assefa M. Melesse Received: 8 May 2016; Accepted: 15 June 2016; Published: 18 June 2016

Abstract: The BeiDou satellite-induced code biases have been confirmed to be orbit type-, frequency-, and elevation-dependent. Such code-phase divergences (code bias variations) severely affect absolute precise applications which use code measurements. To reduce their adverse effects, an improved correction model is proposed in this paper. Different from the model proposed by Wanninger and Beer (2015), more datasets (a time span of almost two years) were used to produce the correction values. More importantly, the stochastic information, i.e., the precision indexes, were given together with correction values in the improved model. However, only correction values were given while the precision indexes were completely missing in the traditional model. With the improved correction model, users may have a better understanding of their corrections, especially the uncertainty of corrections. Thus, it is helpful for refining the stochastic model of code observations. Validation tests in precise point positioning (PPP) reveal that a proper stochastic model is critical. The actual precision of the corrected code observations can be reflected in a more objective manner if the stochastic model of the corrections is taken into account. As a consequence, PPP solutions with the improved model outperforms the traditional one in terms of positioning accuracy, as well as convergence speed. In addition, the Melbourne-Wübbena (MW) combination which serves for ambiguity fixing were verified as well. The uncorrected MW values show strong systematic variations with an amplitude of half a wide-lane cycle, which prevents precise ambiguity determination and successful ambiguity resolution. After application of the code bias correction models, the systematic variations can be greatly removed, and the resulting wide lane ambiguities are more likely to be fixed. Moreover, the code residuals show more reasonable distributions after code bias corrections with either the traditional or the improved model. Keywords: BeiDou; code bias variation; code-phase divergence; stochastic model; precise point positioning

1. Introduction The Chinese BeiDou navigation satellite system (abbreviated as BDS, or BeiDou, or COMPASS) has already launched a regional navigation service by the end of 2012, and continues to develop a global system in the near future. As of October 2015, a total of 20 BeiDou satellites including six geostationary orbit (GEO), seven inclined geosynchronous orbit (IGSO), and seven medium-altitude Earth orbit (MEO) satellites have been launched, transmitting triple-frequency signals centered at B1 (1561.098 MHz), B2 (1207.14 MHz), and B3 (1268.52 MHz). However, it is worth mentioning that only 13 operational BeiDou satellites are presently available (5GEO + 5IGSO + 3MEO) for public use. Sensors 2016, 16, 909; doi:10.3390/s16060909

www.mdpi.com/journal/sensors

Sensors 2016, 16, 909

2 of 16

With the advent of BeiDou, a large amount of research has been dedicated to the aspects of precise orbit and clock determination [1–3], regional ionospheric and tropospheric delay modeling [4–6], differential code bias and timing group delay corrections [7–10], multipath effects analysis [11,12], triple-frequency ambiguity resolution [13,14], and precise positioning [1,15–19], etc. When it comes to the characteristics of BeiDou signals, Nadarajah et al. first discovered the BeiDou Inter-Satellite-Type Biases (ISTBs), and deeply investigated the impacts of such biases on real-time kinematic (RTK) positioning and mixed receiver attitude determination [20,21]. A recent study revealed that the code-phase divergences, which are absent for GPS, GLONASS, and Galileo, however, are commonly found in BeiDou IGSO and MEO satellites [22]. The systematic variation was first mentioned by Hauschild et al. [11], but the origin of the systematic variation was not identified. The same code bias variation was also detected by Montenbruck et al. through the analysis of multipath combination [23,24]. They found that the systematic bias is elevation-dependent, which varies by 0.4–0.6 m from horizon to zenith, particularly pronounced for the B1 signal. Since such systematic bias was consistently observed with other receivers and antennas it was, therefore, attributed to the transmitting satellites. More characteristics of the systematic bias variations were well discussed in Wanninger and Beer [22]. They identified significant differences between two groups of satellites (i.e., IGSO and MEO) and among carrier frequencies (B1, B2, and B3). A correction model was then proposed to accommodate the satellite induced code bias variations, which is definitely an important contribution to the current BeiDou constellation. Nevertheless, the data used for modeling these variations is very limited. Observations from two 10-day intervals in March 2014 (DoY 070–079) and in June 2014 (DoY 160–169) were exclusively used. The long-term stability of the systematic biases could not be identified conclusively from such limited dataset. Consequently, the effectiveness of the mitigating model with fixed correction values may be a debatable issue. In addition, only the correction values on the specific nodes are provided, whereas the stochastic information, i.e., the precision indexes, are completely missing in their correction model. Such precision indexes are, however, important for the stochastic model establishment in precise applications, such as precise point positioning (PPP). Ignoring the stochastic error of the corrections will lead to an overestimation of the actual precision of code observations. This is particularly pronounced in the case of inaccurate corrections are applied. This paper deals with the same systematic bias variations as earlier described by Hauschild et al. [11], Montenbruck et al. [23,24], and Wanninger and Beer [22]. However, we use more datasets with a time span of almost two years (from January 2014 to October 2015) to model such orbit type-, frequency-, and elevation-dependent code biases. Moreover, the stochastic model of the corrections is considered as well. Following an analysis of the code bias variations and a brief introduction of the observation datasets used, we develop an improved correction model, which produce both correction values and their precision indexes. Afterwards, effectiveness validations for the proposed model are presented. Finally, we end with the conclusions. 2. Code Bias Variation Analysis Providing dual- or triple-frequency GNSS observations, multipath (MP) effects are commonly assessed by a linear combination of code and carrier phase observables [12,25]. That is: MPi pi, jq “ Pi ´

f i2 ` f j2 f i2 ´ f j2

λi ϕi `

2 f j2 f i2 ´ f j2

λj ϕj

(1)

where MPi is the multipath combination on frequency i pi, j “ 1, 2, 3, i ‰ jq, P and ϕ represent code range and carrier phase observables, f and λ are, respectively, the frequency and wavelength of the specific signal. If the phase noise and phase multipath are neglected, the MP combination mainly consists of code multipath, a constant ambiguity term which is a combination of the ambiguities of the two phase measurements, a combined signal delay, and thermal noise. Subtraction of the mean value

Sensors 2016, 16, 909

3 of 16

Sensors 2016, 16, 909

3 of 15

from the measurements removes the phase ambiguities, which are constant if there are no cycle slips. The resulting time are series MPslips. is dominated by code andisthermal noise. constant if there no of cycle The resulting timemultipath series of MP dominated by code multipath The MP combination as shown in Equation (1) is an ionosphere-free and geometry-free and thermal noise. combination. The linear coefficients are selected in such a way that ionospheric and tropospheric delay The MP combination as shown in Equation (1) is an ionosphere-free and geometry-free ascombination. well as all geometric contributions (e.g., clocks, orbits, and antenna movements, etc.) cancel out. The linear coefficients are selected in such a way that ionospheric and tropospheric Therefore, canasbeallused with both static and (e.g., kinematic data. However, it is worth mentioning that the delay as it well geometric contributions clocks, orbits, and antenna movements, etc.) cancel absolute values of multipath are unknown. Only variations of code multipath are available from the out. Therefore, it can be used with both static and kinematic data. However, it is worth mentioning MP series. that the absolute values of multipath are unknown. Only variations of code multipath are available Forthe each there are at least two choices of MP combination using different frequencies among from MPcode, series. triple-frequency observations. Take as an of example; the code using multipath on B1 (MP1) can For each code, there are at leastBeiDou two choices MP combination different frequencies triple-frequency Take (B1, BeiDou example; the codecombinations. multipath on B1 (MP1) beamong derived from Equationobservations. (1) using either B2) as or an (B1, B3) frequency However, can be derived from that Equation (1) using either (B1,than B2) or (B1, B3) frequency However, one should be aware the latter one is noisier the former one due combinations. to a larger multiplication one should be aware that the latter one is noisier thanand the former due tofrom a larger factor. Therefore, to keep a lower noise level, both MP1 MP2 areone derived (B1, multiplication B2), and MP3 is factor. Therefore, to keep a lower noise level, both MP1 and MP2 are derived from (B1, B2), and MP3 computed from (B1, B3) combination in this paper. is computed (B1, B3)the combination paper. Shown in from Figure 1 are MP seriesin ofthis BeiDou satellites observed on JFNG (Wuhan, China) with Shown in Figure 1 are the MP series of BeiDou satellites observed on JFNG China) a Trimble Net R9 receiver. BeiDou satellites are grouped by three different orbit (Wuhan, types, i.e., GEO,with IGSO, a Trimble Net R9 receiver. BeiDou satellites are grouped by three different orbit types, i.e., GEO, and MEO. For comparison reason, the MP series of GPS satellites tracked by the same receiver are IGSO, as and MEO. For comparison reason, the MP series of GPS trackedangles by the are same receiver plotted well in Figure 1. As expected, signals observed atsatellites low elevation more likely are plotted as well in Figure 1. As expected, signals observed at low elevation angles are more likely to be affected by multipath error. On the contrary, a high elevation satellite is less likely to cause to be affected by multipath error. On the contrary, a high elevation satellite is less likely to cause multipath error. However, unexpected code biases (drifts), which are absent for GPS, are found in multipath error. However, unexpected code biases (drifts), which are absent for GPS, are found in BeiDou MEO and IGSO satellites. The code bias variation reaches approximately 1 m at high elevation BeiDou MEO and IGSO satellites. The code bias variation reaches approximately 1 m at high elevation angles, particularly for the MEOs on the B1 signal. Similar variation and its characteristics have been angles, particularly for the MEOs on the B1 signal. Similar variation and its characteristics have been well discussed in Wanninger and Beer [22]. Since the code bias variation was consistently observed well discussed in Wanninger and Beer [22]. Since the code bias variation was consistently observed on on other stations equipped with different receivers and antennas, such bias thus was attributed to other stations equipped with different receivers and antennas, such bias thus was attributed to the thetransmitting transmitting satellites. addition, anomalies found MP3 series GPS (i.e.,code code satellites. In In addition, anomalies cancan be be found in in thethe MP3 series of of GPS (i.e., multipath effect on L5 signal), which show significantly larger values at high elevation angles. Similar multipath effect on L5 signal), which show significantly larger values at high elevation angles. Similar code-phase divergences, signalreflections, reflections,had hadbeen been earlier code-phase divergences,which whichwere wereattributed attributed to to satellite-internal satellite-internal signal earlier reported for GPS SVN49 and SVN62 in Spring and Dislssner [26]; Montenbruck et al. [27]; and reported for GPS SVN49 and SVN62 in Spring and Dislssner [26]; Montenbruck et al. [27]; and Hauschild Hauschild al. 2012 [11]. et al. 2012 et [11]. GPS GEO IGSO MEO

2

1 M P2 [m ]

M P 1 [m ]

1 0

0 -1

-1 -2

GPS GEO IGSO MEO

2

0

15

30 45 60 elevation [degrees]

75

90 GPS

2

-2

0

15

30 45 60 elevation [degrees] GEO IGSO MEO

75

90

MP3 [m]

1 0 -1 -2

0

15

30 45 60 elevation [degrees]

75

90

Figure 1. BeiDou and GPS code MP series against elevation angles observed from JFNG (DoY 081/2015). Figure 1. BeiDou and GPS code MP series against elevation angles observed from JFNG (DoY 081/2015).

To investigate the degree of linear dependence between MP and elevation, the Pearson correlation coefficient [28] was applied and computed as:

Sensors 2016, 16, 909

4 of 16

To investigate the degree of linear dependence between MP and elevation, the Pearson correlation coefficient [28] was applied and computed as: ř ρ MP,E “ b

ř

MP ¨ E ´

ř p MP2 ´

MP N

ř

E

ř p MPq2 ř 2 qp E N

´

p

E q2 N q

ř

(2)

where E and N represent satellite elevation angle and the number of samples, respectively. The correlation coefficient ρ MP,E ranges from ´1 to 1. A value of 1 or ´1 implies that a linear equation describes the positive or negative relationship between MP and E perfectly, with all data points lying on a line. A value of 0 implies that there is no linear correlation between the two variables. Table 1 lists the corresponding correlation coefficients for both GPS and BeiDou satellites. In general, the correlation coefficient under 0.3 represents low correlation, while the coefficient between 0.3 and 0.5 indicates normal. If the correlation coefficient exceeding 0.5 but less than 0.7, a medium correlation exists between two components. The correlation coefficient between 0.7 and 0.9 represents a high correlation. Additionally, if the correlation coefficients were beyond 0.9, there would be a strong correlation between two components. The coefficients in Table 1 demonstrate a medium correlation between MP and elevation for BeiDou IGSOs and MEOs, whereas the correlations of GPS and BeiDou GEO satellites are rather weak or even uncorrelated. Nevertheless, one should keep in mind that the GEOs are visible at almost constant elevation angles within restricted region (i.e., Asia-Pacific region). It is unable to infer whether such code biases affect BeiDou GEO code measurements. Hence, the GEOs are excluded in the following context. Table 1. Correlations between MP and elevation series. BeiDou GPS MP1 MP2 MP3

´0.04 0.08 0.12

GEO

IGSO

MEO

´0.01 0.00 ´0.02

´0.49 ´0.48 v0.39

´0.61 ´0.60 ´0.36

3. Observation Datasets As of August 2015, the multi-GNSS experiment (MGEX) [29] offers a global network of approximately 70 stations of GNSS receivers capable to receive BeiDou signals. Most of these stations are equipped with Trimble NetR9 receivers providing triple-frequency BeiDou observations. However, only a few stations are equipped with other receiver types, such as Septentrio POLARX4/POLARX4TR/ASTERX3 and Leica GR25 receivers, providing dual-frequency observations on frequencies B1 and B2. To achieve a global coverage of receiving stations and to be able to compare different receiver types, six stations (CHPG at Cachoeira Paulista, Brazil, CUT0 at Perth, Australia, GMSD at Nakatane town, Japan, JFNG at Wuhan, China, NKLG at Libreville, Gabon, and REUN at Le Tampon, France) equipped with Trimble receivers and four stations (CEBR at Cebreros, Spain, KOUR at Kourou, French Guiana, MAL2 at Malindi, Kenya, and MGUE at Malargue, Argentina) equipped with Septentrio receivers were selected as shown in Figure 2. Data recorded on the first day of each month from January 2014 to October 2015 were used as the core datasets for this study. The datasets spanning almost two years facilitate analysis of a long-term variation of the BeiDou satellite induced code bias. However, it is worth mentioning that some of the data may be unavailable during this period due to station related reasons, e.g., new settled devices or upgrading issues, etc. Whatever the reason, an average of 13–14 days data is usable for each station.

Sensors 2016, 16, 909 Sensors 2016, 16, 909 Sensors 2016, 16, 909

5 of 16 5 of 15 5 of 15

CEBR CEBR KOUR KOUR MGUE MGUE

JFNG JFNG

MAL2 NKLG MAL2 NKLG

GMSD GMSD CUT0 CUT0

CHPG CHPG

REUN REUN

Septentrio Septentrio

Trimble Trimble

Figure 2. Distribution of the selected MGEX stations used in this study. study. Figure 2. Distribution of the selected MGEX stations used in this study.

Figure 333 shows shows the the number number of of MP MP samples samples against against elevation elevation bins bins with with an an interval interval of of 10° 10°˝ for for Figure shows the number of MP samples against elevation bins with an interval of 10 for Figure BeiDou MEO and IGSO satellites. As already mentioned, stations equipped with Septentrio receivers BeiDou MEO MEO and and IGSO IGSO satellites. satellites. As Asalready already mentioned, mentioned, stations stations equipped equipped with with Septentrio Septentrio receivers receivers BeiDou provide only B1 and B2 signals and, thus, the number of MP3 for MEOs is much smaller than that of of provide only B1 and B2 signals and, thus, the number of MP3 for MEOs is much smaller than that of provide only B1 and B2 signals and, thus, the number of MP3 for MEOs is much smaller than that MP1 and MP2. As to IGSOs, the number of MP3 is comparable to that of MP1 and MP2 since the MP1 and of of MP3 is comparable to that of MP1 and MP2 since since the IGSO MP1 and MP2. MP2. As AstotoIGSOs, IGSOs,the thenumber number MP3 is comparable to that of MP1 and MP2 the IGSO satellites are mostly tracked by triple-frequency Trimble receivers in Asia-pacific area. satellites are mostly tracked by triple-frequency Trimble receivers in Asia-pacific area. Furthermore, IGSO satellites are mostly tracked by triple-frequency Trimble receivers in Asia-pacific area. Furthermore, the number ofisMP MP samples is relative relative smallerand at the the lowest and highest highest elevation bins the number ofthe MPnumber samples relative smaller at the smaller lowest highest elevation binselevation due to spatial Furthermore, of samples is at lowest and bins due to spatial configuration reason. configuration reason. due to spatial configuration reason.

Num Number berofofs am s amples ples

6 6 4 4

4

B1 B1 B2 B2 B3 B3

2 2 0 0

5 5

x 104 15 x 10 15

MEO MEO

15 25 35 45 55 65 75 85 15 25 35 45 55 65 75 85 elevation [degrees] elevation [degrees]

Num Number berofofs am s amples ples

4

x 104 8 x 10 8

IGSO IGSO B1 B1 B2 B2 B3 B3

10 10 5 5 0 0

5 5

15 25 35 45 55 65 75 85 15 25 35 45 55 65 75 85 elevation [degrees] elevation [degrees]

Figure 3. Number of MP samples against elevation bins with an interval of 10° for BeiDou MEO (left) Figure Figure 3. 3. Number Number of of MP MP samples samples against against elevation elevation bins bins with with an an interval interval of of 10° 10˝ for forBeiDou BeiDouMEO MEO(left) (left) and IGSO satellites (right). and IGSO satellites (right). and IGSO satellites (right).

4. Improved Correction Correction Model 4. 4. Improved Improved CorrectionModel Model As shown shown in in Figure Figure 1, 1, the the code code bias bias varies varies linearly linearly with with elevation. elevation. To To obtain obtain the the best best fitting fitting As As shown in Figure 1, the code bias varies linearly with elevation. To obtain the best fitting results, results, continuous continuous piecewise piecewise linear linear functions functions were were used used to to model model the the code-phase code-phase divergences. divergences. Note Note results, continuous piecewise linear functions were used to model the code-phase divergences. Noteshare that the that the term continuous is used in the sense that the adjacent segments of the function the that term continuous is the used in the that the segments adjacent segments of theshare function share end the termthe continuous is used in sense thatsense the adjacent of the function the same same end point. point. Given aa set set of of MP MP series, series, ((E which can can be be divided divided into into several several consecutive consecutive Ei,,MP MPi ) ,, which same point.end Given a setGiven of MP series, pEi , MPi q, which divided into several consecutive segments, the i cani )be segments, objective the following following objective functions should be satisfied satisfied tofits: find the the best best fits: fits: following functions should be satisfied to be find the bestto segments, the objective functions should find m 1 nn j 1 1 jn j m m´

ř ř ( f  MP )22  min  j ,i) MPjj,i min SSS “   ( pffjjj,i,,ii ´MP  min ,i q “ j“ 11 ii“ 1 jj 1 i 1 1

q ´ f pa q “ 0 s.t. s.t.ff j((afaj pa j )j f j j1`(1a j )j  0 s.t. j j )  f j 1 ( a j )  0 where where where ‚ ii“ 1, 2, ¨ ¨ ¨ n: point index; 1,2, nnumber point index; index; n:n numberofof ofpoints points(samples); (samples); i  1,2, nn:: point :: number points (samples);  ‚ j “ 1, 2, ¨ ¨ ¨ m ´ 1: segment index; m: number of nodes; j 1,2, 1,2, m m11:: segment m:: number segment index; index; m number of of nodes; ‚ ajj :elevation angle on a specific node (i.e., the end nodes; point of j-th segment); elevation angle angle on on aa specific specific node node (i.e., (i.e., the the end end point point of of j-th j-th segment); segment); aajj:: elevation 

(3) (3) (3)

Sensors 2016, 16, 909 Sensors 2016, 16, 909 Sensors 2016, 16, 909

6 of 16 6 of 15 6 of 15

bj :j MP : MPvalue valueonona aspecific specificnode; node; bb j : MP value on a specific node; f j,i function used in the j-th j-th segment. m q: f “f j pE f (j,iE, b1 , bb2 , ¨b¨ ¨,b b fitted ) : fitted function used in the segment. f jj,,ii  f jj ( E jj,,ii , b11, b22 ,bmm ) : fitted function used in the j-th segment. Solving Solving the the normal normal linear linear system system of of Equation Equation (3) (3) with with least least squares squares method method produces produces the the MP MP Solving the normal linear system of Equation (3) with least squares method produces the MP values values on on all all nodes. nodes. Then, Then, the the precision precision of of estimates, estimates, i.e., i.e., the the root root mean mean square square(RMS) (RMS)isisderived derivedas: as: values on all nodes. Then, the precision of estimates, i.e., the root mean square (RMS) is derived as: g f nnj j f nř f j (p ff j,i ´  MP MPj,ij ,iq)222 f ( f j ,i  MP ) (4) e i “1 (4) RMS jj “ RMS  i 11 nj ,i ´ 1 j ,i (4) RMS  nj  1 ‚  ‚ 

 

j

n jj  1

To To ensure ensure enough enough samples samples are are used used and and to to get get reliable reliable fitting fitting functions, functions, we we selected selected the the nodes nodes ˝ To ensure enough samples are used and to get reliable fitting functions, we selected the nodes to be separated separatedbyby10° 10of elevation, of elevation, the node the elevation angles range to be and and the node valuesvalues for thefor elevation angles range from 5° tofrom 85°. ˝ be ˝ . Representative separated by 10° of elevation, and the node values for the elevation angles range from 5° to 85°. 5to to 85 results of such continuous piecewise linear models are presented Representative results of such continuous piecewise linear models are presented in Figures 4 and in 5. Representative of4 such continuous piecewise linear modelsstation are presented in 4 China, and 5. Figures 5. results Figure shows the results forstation JFNG, a MGEX in Figures Wuhan, Figure 44 and shows the results for JFNG, a MGEX located in Wuhan, located China, which is equipped Figure is4 equipped shows thewith results for JFNG, a MGEX station located in Wuhan,data China, which istest equipped which a Trimble receiver. Alldata available during period with a Trimble NetR9 receiver. AllNetR9 available BeiDou duringBeiDou the test period (i.e.,the January 2014– withJanuary a Trimble NetR9 receiver. All available BeiDou data during the test period (i.e.,B3 January 2014– (i.e., 2014–October 2015) were used to model the code biases on B1, B2, and frequencies. October 2015) were used to model the code biases on B1, B2, and B3 frequencies. Scatterplots on each October 2015) used to model the biasesof ondifferent B1, B2, and B3which frequencies. Scatterplots on each Scatterplots onwere each node represent the code MP days, values reflects the bias variation node represent the MP values of different which clearlydays, reflects theclearly variation of code with node represent the MP values of different days, which clearly reflects the variation of code bias with of code bias with time. The solid lines connect the averaged MP values on consecutive segments. time. The solid lines connect the averaged MP values on consecutive segments. Likewise, the code time. The the solidcode linesbiases connect the averaged MP on stations consecutive segments. Likewise, the code Likewise, from all thevalues selected the same day (DoY biases estimated from allestimated the selected stations on the same dayon (DoY 001/2015) were 001/2015) shown in biases estimated from all the selected stations on the same day (DoY 001/2015) were shown in were shown in Figurescatter 5. Differently, plots on represent the same node represent MP values of different Figure 5. Differently, plots onscatter the same node MP values of different stations, which Figure 5. Differently, scatter plots on the same node represent MP values of different stations, which stations, gives impression us an intuitive impression of their consistency. gives us which an intuitive of their consistency. gives us an intuitive impression of their consistency. MEO MEO

B1 B2 B3 B1 B2 B3

0 0 -0.5 -0.5 -1 -1

0.4 0.4

IGSO IGSO

0.2 0.2 MP [m] MP [m]

MP [m] MP [m]

0.5 0.5

B1 B2 B3 B1 B2 B3

0 0 -0.2 -0.2

-0.4 15 25 35 45 55 65 75 85 -0.4 5 15 25 35 45 55 65 75 85 15 25 35 45 55 65 75 85 5 15 25 35 45 55 65 75 85 elevation [degrees] elevation [degrees] elevation [degrees] elevation [degrees] Figure 4. Piecewise Piecewise linear MP models for BeiDou MEO MEO (left) (left) and and IGSO satellites satellites (right) on on JFNG Figure linear MP MP models models for BeiDou BeiDou Figure 4. 4. Piecewise linear for MEO (left) and IGSO IGSO satellites (right) (right) on JFNG JFNG station using data from different days (during January 2014–October 2015). station station using using data data from from different different days days (during (during January January 2014–October 2014–October2015). 2015). 5 5

MEO MEO

0 0 -0.5 -0.5 -1 -1

0.4 0.4

B1 B2 B3 B1 B2 B3

0.2 0.2 MP [m] MP [m]

MP [m] MP [m]

0.5 0.5

IGSO IGSO

B1 B2 B3 B1 B2 B3

0 0 -0.2 -0.2

-0.4 15 25 35 45 55 65 75 85 -0.4 5 15 25 35 45 55 65 75 85 15 25 35 45 55 65 75 85 5 15 25 35 45 55 65 75 85 elevation [degrees] elevation [degrees] elevation [degrees] elevation [degrees] Figure 5. Piecewise linear MP models for BeiDou MEO (left) and IGSO satellites (right) using data Figure 5. 5. Piecewise Piecewise linear linear MP MP models models for for BeiDou MEO MEO (left) (left) and and IGSO IGSO satellites satellites (right) (right) using using data data Figure from different stations (equipped with six BeiDou Trimble receivers and four Septentrio receivers). from different different stations stations (equipped (equipped with with six six Trimble Trimblereceivers receiversand andfour fourSeptentrio Septentrioreceivers). receivers). from 5 5

Obviously, such code biases are orbit type-, elevation-, and frequency-dependent, which is Obviously, such code biases are orbit type-, elevation-, and frequency-dependent, which is consistent with Wanninger and Beer [22]. The code measurements of MEO satellites tracked on B1 consistent with Wanninger and Beer [22]. The code measurements of MEO satellites tracked on B1 frequency at high elevations are seriously affected by the code bias, which reaches over 0.5 m when frequency at high elevations are seriously affected by the code bias, which reaches over 0.5 m when

Sensors 2016, 16, 909

7 of 16

1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 -0.2 -0.2 -0.4 -0.4

B1 B1 B2 B2 B3 B3

5 5

15 15

MEO MEO

Correction values [m][m] Correction values

Correction values [m][m] Correction values

Obviously, such code biases are orbit type-, elevation-, and frequency-dependent, which is consistent with Wanninger and Beer [22]. The code measurements of MEO satellites tracked on B1 Sensors 2016, 16, 909 7 of 15 frequency at high Sensors 2016, 16, 909 elevations are seriously affected by the code bias, which reaches over 0.5 m when 7 of 15 ˝ The code biases on B2 and B3, however, are much smaller than those on the the elevation elevation rises rises up up to to 70 70°.. The code biases on B2 and B3, however, are much smaller than those on the elevation rises up to 70°. The code similar biases on B2 and however, are much smaller thanHowever, those on B1. As to to the the IGSO IGSO group of satellites, satellites, effects canB3, befound found among triplefrequencies. frequencies. B1. As group of similar effects can be among triple However, B1. As to the IGSO group of satellites, similar effects can be found among triple frequencies. However, the are much smaller compared to the values values of of code codebias biasas aswell wellas astheir theirfrequency-dependent frequency-dependentdifferences differences are much smaller compared the values of code bias as well as their frequency-dependent differences are much smaller compared those of MEOs. In general, the code biases show a good consistency among different stations (equipped to those of MEOs. In general, the code biases show a good consistency among different stations to those of MEOs. In general, the code biases a good consistency among different stations with different receivers) different days and,show thus, make possible to be correctly modeled with (equipped with differentand receivers) and different days and,it thus, make it possible to be correctly (equipped with different receivers) and different days differences and, thus, do make it in possible to beif correctly limited data. Nevertheless, it is worth noting that slight exist the models different modeled with limited data. Nevertheless, it is worth noting that slight differences do exist in the modeled with limited data. Nevertheless, it is worth noting that slight differences do exist in the samples used. models ifare different samples are used. models if different samples are used. best fitting model of code biases, MP series derived from all the Therefore, in order Therefore, in order to to obtain obtain the the best fitting model of code biases, MP series derived from all the Therefore, induring order to obtain of the best fitting modelwere of code biases, MP seriescorrection derived from all the selected used selected stations stations during aa period period of almost almost two two years years were used to to produce produce the the correction model model for for selected stationsinduced during acode period of In almost two to years were usedvalues, to produce the correction precisions model for BeiDou satellite bias. addition the correction the corresponding BeiDou satellite induced code bias. In addition to the correction values, the corresponding precisions BeiDou code bias. In addition to the correction the corresponding precisions are essential for to aa better of their Notice there are also also satellite essentialinduced for the the users users to have have better understanding understanding ofvalues, their corrections. corrections. Notice that that there is is are also essential for theatusers to have afor better understanding of their corrections. Notice that thereare is very large RMS-values low elevation the reason that signals observed at low elevation angles very large RMS-values at low elevation for the reason that signals observed at low elevation angles very large RMS-values atby low elevationerror for the reason thatwhich signals observed atprecision low elevation angles more likelylikely to be to affected multipath stations, decrease estimation. are more be affected by multipathfrom error from stations, which the decrease theofprecision of ˝ are more likely to bethe affected bycorrection multipathvalues error and fromtheir stations, which decrease thefrom precision of˝ Figures 6 and 7 show resulting RMS for elevations range 5 to 85 estimation. Figures 6 and 7 show the resulting correction values and their RMS for elevations range ˝ . Thethe estimation. Figures 6 and show resulting correction values and their RMS for elevations range with separation of 710 final parameters areparameters summarized Table 2. from a5°node to 85° with a node separation ofmodel 10°. The final model areinsummarized in Table 2. from 5° to 85° with a node separation of 10°. The final model parameters are summarized in Table 2.

25 35 45 55 65 25elevation 35 45 55 65 [degrees] elevation [degrees]

75 75

85 85

1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 -0.2 -0.2 -0.4 -0.4

IGSO IGSO B1 B1 B2 B2 B3 B3

5 5

15 15

25 35 45 55 65 25elevation 35 45 55 65 [degrees] elevation [degrees]

75 75

85 85

MEO MEO

0.8 0.8

B1 B1 B2 B2 B3 B3

0.6 0.6 0.4 0.4 0.2 0.2

RMRM SS of of corrections [m[m ] ] corrections

RMRM SS of of corrections [m[m ] ] corrections

Figure 6. Suggestedcorrection correction values of of code code bias bias (orbit (orbit type-, frequency-, and elevation-dependent). Figure type-, frequency-, frequency-,and andelevation-dependent). elevation-dependent). Figure6.6.Suggested Suggested correction values values of code bias (orbit type-, IGSO IGSO

0.8 0.8

B1 B1 B2 B2 B3 B3

0.6 0.6 0.4 0.4 0.2 0.2

0 0 5 15 25 35 45 55 65 75 85 15 25 35 45 55 65 75 85 5 15 25elevation 35 45[degrees] 55 65 75 85 15 25elevation 35 45[degrees] 55 65 75 85 elevation [degrees] elevation [degrees] Figure 7. Suggested RMS of correction values (orbit type-, frequency-, and elevation-dependent). Figure Figure 7. 7. Suggested Suggested RMS RMS of of correction correction values (orbit type-, frequency-, and elevation-dependent).

0 0

5 5

Table 2. Correction values and their precisions for BeiDou MEO/IGSO code measurements. Table 2. Correction values and their precisions for BeiDou MEO/IGSO code measurements.

To make use of them, correction values for elevations betweenRMS 5˝ and 85˝ can be obtained by Correction Values (m) of Corrections (m) Elevation Correction Values (m)IGSOs RMS of corrections Corrections (m) linear interpolation from the nearest two nodes (Equation (5)), whereas the on nodes of 5˝ MEOs MEOs IGSOs Elevation Nodes (°) ˝ MEOs IGSOs MEOs IGSOs B1 used B2for the B3elevation B1 binsB2 B1 B2 B3 the MP B1 correction B2 B3 andNodes 85 are of (0,5)B3 and (85,90). Assuming values (°) simply B1 B2 B3 B1 B2 B3 B1 B2 B3 B1 B2 B3 −0.140 −0.060 −0.148 0.721 0.588 0.580 0.709 at node5 points−0.109 are uncorrelated and−0.101 neglecting the −0.065 interpolation error, the precision of 0.564 a given0.576 point 5 −0.109 0.588 0.580 0.709 0.564 0.576 15 −0.169 −0.140 −0.148 −0.060 −0.087 −0.101 −0.203 −0.148 −0.250 −0.065 −0.162 0.721 0.605 0.480 0.499 0.651 0.532 0.582 can be15determined the law of propagation of variance Equation −0.169 by−0.148 −0.087 −0.203 −0.250 −0.162 in0.605 0.480(5):0.499 0.651 0.532 0.582 25 25 35 35 45 45 55 55 65 65 75 75 85 85

−0.150 −0.150 −0.105 −0.105 0.004 0.004 0.181 0.181 0.411 0.411 0.674 0.674 0.853 0.853

−0.121 −0.121 −0.062 −0.062 0.047 0.047 0.185 0.185 0.326 0.326 0.477 0.477 0.600 0.600

−0.070 −0.070 −0.053 −0.053 0.022 0.022 0.096 0.096 0.180 0.180 0.280 0.280 0.373 0.373

−0.222 −0.222 −0.123 −0.123 −0.066 −0.066 0.036 0.036 0.107 0.107 0.163 0.163 0.245 0.245

−0.224 −0.224 −0.110 −0.110 −0.043 −0.043 0.044 0.044 0.106 0.106 0.178 0.178 0.260 0.260

−0.168 −0.168 −0.078 −0.078 −0.049 −0.049 0.021 0.021 0.068 0.068 0.130 0.130 0.208 0.208

0.476 0.476 0.388 0.388 0.333 0.333 0.293 0.293 0.275 0.275 0.261 0.261 0.233 0.233

0.373 0.373 0.291 0.291 0.254 0.254 0.220 0.220 0.194 0.194 0.188 0.188 0.173 0.173

0.401 0.401 0.290 0.290 0.258 0.258 0.241 0.241 0.211 0.211 0.206 0.206 0.198 0.198

0.500 0.500 0.403 0.403 0.389 0.389 0.308 0.308 0.262 0.262 0.251 0.251 0.217 0.217

0.371 0.371 0.297 0.297 0.278 0.278 0.230 0.230 0.210 0.210 0.213 0.213 0.195 0.195

0.409 0.409 0.303 0.303 0.244 0.244 0.223 0.223 0.208 0.208 0.212 0.212 0.190 0.190

Sensors 2016, 16, 909

8 of 16

MPpeq “ MP0 ` pMP1 ´ MP0 q ¨ Ee´´EE00 1

2

(5)

2

2 σMP “ p EE1´´Ee0 q ¨ σ02 ` p Ee´´EE00 q ¨ σ12 peq 1

1

where Ej , MP and σj2 pj “ 1, 2q are, respectively, the elevation angles, MP values and the corresponding Sensors 2016, 16, j909 8 of 15 RMS on the nearest two nodes; MPpeq and σMPpeq are the interpolated MP correction and RMS at the givenTo elevation angle e. make use of them, correction values for elevations between 5° and 85° can be obtained by linear interpolation from the nearest two nodes (Equation (5)), whereas the corrections on nodes of Table Correction andelevation their precisions for(0,5) BeiDou code measurements. 5° and 85° are 2.simply usedvalues for the bins of andMEO/IGSO (85,90). Assuming the MP correction values at node points are uncorrelated and neglecting the interpolation error, the precision of a given Correction Values (m) RMS of Corrections (m) point can be determined by the law of propagation of variance in Equation (5): Elevation MEOs

Nodes (˝ ) B1

IGSOs

MEOs

e B1  E0 B2 B1 B2 B3 MPB3(e)  MP 0  ( MP1  MP0 )  ´0.060 ´0.101 ´0.148 ´0.065 E0.721 1  E0 0.588

B2

IGSOs

B3

B1

B2

B3

5 ´0.109 ´0.140 0.580 0.709 0.564 0.576 (5) 15 ´0.169 ´0.148 ´0.087 ´0.203 ´0.250 ´0.162 0.605 0.480 0.499 0.651 0.532 0.582 e  E0.476 E1  e´0.224 2 2 2 ´0.168 2 0 2 0.373 25 ´0.150 ´0.121 ´0.070 ´0.222 0.401 0.500 0.371 0.409 )  0  ( ) 1 MP ( e )  ( 35 ´0.105 ´0.062 ´0.053 ´0.123 E1  E0´0.110 ´0.078 E1  0.388 E0 0.291 0.290 0.403 0.297 0.303 45 0.004 0.047 0.022 ´0.066 ´0.043 ´0.049 0.333 0.254 0.258 0.389 0.278 0.244 55 0.181 0.1852 0.096 0.036 0.044 0.021 0.293 0.220 0.241 0.308 0.230 0.223  j ( j0.180  1, 2) 0.107 are, respectively, elevation values the where65 E j , MP j and0.326 0.411 0.106 0.068the 0.275 0.194angles, 0.211 MP 0.262 0.210and 0.208 75 0.674 0.477 0.280 0.163 0.178 0.130 0.261 0.188 0.206 0.251 0.213 0.212 (e) and corresponding RMS on the nearest two interpolated MP 85 0.853 0.600 0.373 0.245 nodes; 0.260 MP0.208 0.233 MP 0.173 0.198the0.217 0.195 0.190 ( e ) are

Differences [m]

correction and RMS at the given elevation angle e . Figure 8 8 shows shows the correction differences and Beer Beer [22]. [22]. In Figure the correction differences with with respect respect to to Wanninger Wanninger and In general, general, the differences are mostly within 0.1 m, which implies that our model agrees well with the existing the differences are mostly within 0.1 m, which implies that our model agrees well with the existing one. one. However, the values of the differences at low elevation bins are much bigger, which reach However, the values of the differences at low elevation bins are much bigger, which reach 0.2–0.4 m. 0.2–0.4 This by is caused by the significantly larger noise levelelevation at low elevation The precision This is m. caused the significantly larger noise level at low angles. angles. The precision of MP of MP corrections asin listed 2 confirms corrections as listed Tablein2Table confirms this. this. 0.4 0.3 0.2 0.1

B1(IGSO) B2(IGSO) B3(IGSO) B1(MEO) B2(MEO) B3(MEO)

0 -0.1 -0.2 -0.3 -0.4

5

15

25 35 45 55 65 elevation [degrees]

75

85

Figure8.8.Differences Differences of of the the code bias correction values Figure values compared comparedwith withthose thoseof ofWanninger Wanningerand andBeer. Beer.

5. Validation 5. Validationof of the the Models Models To verify verifythe theproposed proposed correction model, absolute positioning techniques, as standard To correction model, absolute positioning techniques, such such as standard point point positioning (SPP) and PPP, are good choices since they both are affected by the satellite-induced positioning (SPP) and PPP, are good choices since they both are affected by the satellite-induced code biases code biases may be be obscured by code biases described described above. above.ItItshould shouldbebenoted notedthat, that,however, however,such such code biases may obscured the residual ionospheric delay and broadcast ephemeris (i.e., satellite orbits and clock corrections) errors, by the residual ionospheric delay and broadcast ephemeris (i.e., satellite orbits and clock corrections) which which are regarded as theasdominating errors in SPP. Consequently, the errors, are regarded the dominating errors in SPP. Consequently, thecode codebias bias correction correction has has marginal effect on the the positioning positioning results. results. To clearly demonstrate demonstrate the the effectiveness effectiveness of of the the improved improved marginal effect on To clearly correction model, were performed withwith three different processing schemes. In the first scheme, correction model,PPP PPPtests tests were performed three different processing schemes. In the first the satellite-induced code biases were kept without any treatment. In the second scheme, the traditional scheme, the satellite-induced code biases were kept without any treatment. In the second scheme, correction model provided in Wanninger and Beer [22] and wereBeer used to were remove the biases. The the traditional correction model provided in Wanninger [22] used to code remove the code improved correction model which contains correction values as well as precision indexes were applied for mitigating the same code biases in the last scheme. In the interest of brevity, “None”, “Traditional”, and “Improved” are used to denote the above-mentioned three different schemes throughout the rest of this article, if there is no additional explanation. The common processing strategies for PPP are summarized in Table 3. Precise BeiDou orbits and clock corrections provided

Sensors 2016, 16, 909

9 of 16

biases. The improved correction model which contains correction values as well as precision indexes were applied for mitigating the same code biases in the last scheme. In the interest of brevity, “None”, “Traditional”, and “Improved” are used to denote the above-mentioned three different schemes throughout the rest of this article, if there is no additional explanation. The common processing strategies for PPP are summarized in Table 3. Precise BeiDou orbits and clock corrections provided by German Research Center for Geoscience (GFZ, Potsdam, Germany) were used. The positioning results were compared with the IGS weekly solutions [30]. In general, the reference coordinates have an accuracy of few millimeters. Table 3. PPP processing strategy. Items

Models

Estimator (engine)

Kalman filter, TriP software [31]

Observations

B1/B2 code and carrier phase measurements

Sampling rate

30 s

Elevation cutoff

10˝

Weighting scheme

Elevation dependent; 3 mm and 0.3 m for raw phase and code, respectively

Ionospheric delay

Eliminated by Ionosphere-free combination (s)

Tropospheric delay

Dry component: corrected with GPT model [32] wet component: estimated as random-walk process, GMF mapping function

Relativistic Effect

Applied

Station displacement

Corrected by IERS Convention 2010, including Solid Earth tide and ocean tide loading [33]

Satellite antenna phase center offset

Corrected with conventional PCO values from MGEX [34]

Receiver antenna phase center offset

Corrected by the same PCO values as GPS

Phase-windup effect

Corrected [35]

Receiver clock

Estimated, epoch-wise solution

Station coordinate

Estimated, epoch-wise solution

Phase ambiguities

Estimated, constant for each arc; float value

5.1. Positioning Error The validation of the correction is verified by PPP solution using five stations: CUT0, JFNG, XMIS (Christmas Island, Australia), GMSD (Nakatane town, Japan), NNOR (New Norcia, Australia), during days 45–50 in 2014 and days 130–135 in 2015 while the XMIS and SIN1 (Singapore) stations were excluded from the code bias modeling. Two representative PPP solutions on different stations and different days are shown in Figures 9 and 10. Figure 9 shows the epoch-wise positioning error of CUT0 and JFNG recorded on DoY 047/2014, and Figure 10 shows the positioning error of CUT0 and XMIS recorded on DoY 130/2015. Owing to the highly-weighted carrier phase, the positioning results are not seriously affected by the un-modeled code biases. In general, an accuracy of 0.1 in horizontal and 0.2 m vertically can be achieved after convergence with the current BeiDou constellation even if we leave the code biases uncorrected. By comparing the first two schemes, the code bias corrections with the traditional model does not really yield noticeable improvements in terms of positioning. Sometimes inaccurate correction values may even degrade the positioning accuracy. However, the third scheme outperforms the former two in terms of positioning accuracy as well as convergence speed. This implies that the improved model, which provides not only correction values but also stochastic information, greatly enhances PPP performance, particularly at the initial stage.

Sensors 2016, 16, 2016, 909 16, 909 Sensors Sensors 2016, 16, 909

UpUp [m][m]

North [m][m] North

East [m][m] East

10 of 16

0.2 0.2 0 0 -0.2 -0.2 0.2 0.2 0 0 -0.2 -0.2 0.4 0.4 0 0 -0.4 -0.4

CUT0 CUT0

0 4 0 4

8 12 16 8UTC 12[h] 16 UTC [h] None None

JFNG JFNG

20 24 0 4 20 24 0 4 Traditional Traditional

8 12 16 20 24 8UTC 12[h] 16 20 24 UTC [h] Improved Improved

UpUp [m][m]

North [m][m] North

East [m][m] East

Figure 9. Epoch-wise PPP solutions of CUT0 and JFNG (DoY 047/2014). Figure 9. Epoch-wise PPPsolutions solutions of and JFNG (DoY(DoY 047/2014). Figure 9. Epoch-wise PPP ofCUT0 CUT0 and JFNG 047/2014). 0.2 0.2 0 0 -0.2 -0.2 0.2 0.2 0 0 -0.2 -0.2 0.4 0.4 0 0 -0.4 -0.4

CUT0 CUT0

0 4 0 4

8 12 16 8UTC 12[h] 16 UTC [h] None None

XMIS XMIS

20 24 0 4 20 24 0 4 Traditional Traditional

8 12 16 20 24 8UTC 12[h] 16 20 24 UTC [h] Improved Improved

Figure 10. Epoch-wise PPPsolutions solutions of and XMIS (DoY(DoY 130/2015). Figure 10. Epoch-wise PPP ofCUT0 CUT0 and XMIS 130/2015). Figure 10. Epoch-wise PPP solutions of CUT0 and XMIS (DoY 130/2015).

3D3D position error [m] position error [m]

3D3Dpos ition [m[m ] ] pos itionerror error

JFNG To gain more insightCUT0 into the PPP convergence, the three dimensional (3D) positioning errors of 2 2 JFNG CUT0 the first few hours on CUT0 and JFNG are shown in Figure 11. At the initial stage, large differences 2 2 None None can be observed for both CUT0 and JFNG, positioning results are quiteNone sensitive to None indicating that the 1.5 1.5 Traditional Traditional 1.5 the 1.5 different processing schemes. After a short time convergence (15–30 min), the positioning accuracy Traditional Traditional Improved Improved of CUT0 code bias corrections 1with the traditional model, Improved whereas the Improved 1 is noticeably improved by the 1 time means obtaining a 3D positioning 1 is totally different for JFNG. Note that, the convergence situation error and Zhang [17]. The convergence time for 0.5less than one decimeter, which has been adopted by Li0.5 0.5 0.5 five stations under three schemes and the average convergence time are shown in Table 4. In 0 fact, the correction values may be not accurate enough. 0 The inaccurate corrections may even 00:00 00:30 of code 01:00 01:30 However, 02:00 as already 0 the precision 0 mentioned, 00:00 00:30 the01:00 02:00 degrade observations. stochastic01:30 information 00:00 00:30 01:00 01:30 02:00 00:00 00:30 UTC 01:00 01:30 02:00 UTC [h] [h] of the code bias corrections are completely missing in the traditional model. In other words, the UTC [h] UTC [h] precision of the code observables has been overestimated in this case. This will, unfortunately, lead Figure 11. 3D positioning error of the first two hours on CUT0 and JFNG (DoY 047/2014). to Figure an uncertainty of their corrections, andfirst eventually the on positioning accuracy cannot be always 11. 3D positioning error of the two hours CUT0 and JFNG (DoY 047/2014). enhanced. Compared against the former two schemes, solutions with the improved model show the

Table 4. Convergence time under three schemes for five stations. Table 4. Convergence time under three schemes for five stations.

Convergence Time (min) Convergence Time (min) GMSD

Uncorrected Uncorrected 125

Schemes Schemes Traditional Traditional 108.5

Improved Improved 29.5

10 of 15 10 of 15

North [m]

-0.2 0.2 0 -0.2

Sensors 2016, 16, 909

11 of 16 Up [m]

0.4 0

best positioning performance with an accuracy of ~0.1 m, suggesting that a proper stochastic model -0.4 for the corrections is critical. Reviewing the RMS of the code bias corrections in Table 2, we can find 0 4 8 12 16 20 24 0 4 8 12 16 20 24 that the precision of the provided corrections m. [h] Such an uncertainty is comparable to UTC [h]range 0.2–0.7UTC the assumed noise level of code measurements. Therefore, the actual None Traditional Improvedprecision of the corrected code observations can be reflected in an objective manner once the stochastic model of the corrections is taken into account. Figure 10. Epoch-wise PPP solutions of CUT0 and XMIS (DoY 130/2015). CUT0 None Traditional Improved

1.5 1 0.5 0 00:00

00:30

01:00 UTC [h]

JFNG

2

01:30

3D pos ition error [m ]

3D position error [m]

2

None Traditional Improved

1.5 1 0.5

02:00

0 00:00

00:30

01:00 UTC [h]

01:30

02:00

Figure 11. 11. 3D 3D positioning twohours hourson onCUT0 CUT0 and JFNG (DoY 047/2014). Figure positioningerror errorofofthe the first first two and JFNG (DoY 047/2014). Table 4. Convergence threeschemes schemesforfor five stations. Table 4. Convergencetime time under under three five stations.

Convergence Time (min) Convergence Time (min) GMSD GMSD CUT0 CUT0 JFNG JFNG XIMS XIMS NNOR NNOR Average Average

Schemes Schemes Uncorrected Uncorrected 125 125 130 130 7676 157 157 130 130 123.6 123.6

Traditional Improved Improved Traditional 108.5

108.5 134 134 136 136 162 162 127 127 133.5

133.5

29.5

29.5 3636 32.5 32.5 55 55 33.5 33.5 37.3

37.3

5.2. Wide-Lane Ambiguity In addition to the positioning error, the wide-lane ambiguities, i.e., the so-called Melbourne-Wübbena (MW) linear combination [36,37], were derived from the B1/B2 dual-frequency code and carrier phase measurements. Figures 12 and 13 shows the WL ambiguity and EWL ambiguity of representative IGSO and MEO satellites with different code bias correction schemes. As shown in the upper two plots of Figures 12 and 13, the WL and EWL values derived from original code observables reveal strong elevation-dependent systematic variations with amplitudes of half a wide-lane cycle, which prevents precise ambiguity determination and successful ambiguity resolution. After application of the code bias correction models, the systematic variations have been greatly removed, and the resulting MW series run much more stable and closer to the nearest integers. The fixing rate of MW ambiguities increases from 80.4% to 91.8%. Additionally, the MW values of all involved IGSO (C6, C7, C8, C9, and C10) and MEO (C11, C12, and C14) satellites were calculated and the statistics were given in Figure 14. Likewise, the MW derivations from the latter two schemes agree well with each other. This is reasonable since they have consistent code bias correction values as shown in Figure 8. Moreover, the distributions of the corrected MW values are much more concentrated compared to the one calculated from uncorrected code measurements. Obviously, the corrected MW values are more likely to be fixed, which definitely contribute to the subsequent undifferenced ambiguity resolution.

given in Figure 14. Likewise, the MW derivations from the latter two schemes agree well with each other. This is reasonable since they have consistent code bias correction values as shown in Figure 8. Moreover, the distributions of the corrected MW values are much more concentrated compared to the one calculated from uncorrected code measurements. Obviously, the corrected MW values are more likely to be fixed, which definitely contribute to the subsequent undifferenced ambiguity Sensors 2016, 16, 909 12 of 16 resolution.

MW [cyc]

MW [cyc]

90

1

60

0 -1

30

None

0

1 0 -1

Traditional

1 0 -1 0

Improved 4

8 12 16 UTC [h]

18

20 22 UTC [h]

24 12 of 15 12 of 15

Figure Figure12. 12.Time Timeseries seriesof ofthe thederived derivedMW MWvalue valueon onXMIS XMIS(DoY (DoY130/2015). 130/2015). C11(MEO) C11(MEO)

EML(CYC) EML(CYC)

C09(IGSO) C09(IGSO)

1 1

60 60

0 0 -1 -1

90 90

30 30

NONE NONE

0 0

Elevation(deg) Elevation(deg)

Sensors 2016, 16, 909 Sensors 2016, 16, 909

C11(MEO) Elevation [deg]

MW [cyc]

C9(IGSO)

1 1 0 0 -1 -1

Traditional Traditional

1 1 0 0 -1 -1 0 0

Improved Improved

4 4

8 8

12 12

UTC(h) UTC(h)

16 16

18 18

20 20

22 22

UTC(h) UTC(h)

24 24

Figure 13. 13. Time Time series series of of the thederived derived EWL EWL value valueon onCUT0 CUT0(DoY (DoY122/2015). 122/2015). Figure 13. Time series of the derived EWL value Figure on CUT0 (DoY 122/2015). None None

Traditional Traditional

Improved Improved

10 10 5 5 0 0

15 15 Percentage Percentage[%] [%]

Percentage Percentage[%] [%]

15 15

-0.5 0 0.5 -0.5 0 0.5

-0.5 0 0.5 -0.5 0 0.5 MW [cyc] MW [cyc]

-0.5 0 0.5 -0.5 0 0.5

None None

Traditional Traditional

Improved Improved

-0.5 0 0.5 -0.5 0 0.5

-0.5 0 0.5 -0.5 0 0.5 MW [cyc] MW [cyc]

-0.5 0 0.5 -0.5 0 0.5

10 10 5 5 0 0

Figure14. 14. Histograms ofall allof the derived MWvalues values (left: IGSOs; right: MEOs)right: onXMIS XMIS (DoY130/2015). 130/2015). 14.Histograms Histograms all the derived MW(left: values (left: IGSOs; MEOs) on XMIS Figure of the derived MW IGSOs; right: MEOs) on (DoY (DoY 130/2015).

5.3. Code Code Residual Residual 5.3. Sincethe theun-modeled un-modeledcode codebiases biaseswill willbe bepartly partlyreflected reflectedin inthe theresiduals, residuals,the thecode coderesiduals residuals(i.e., (i.e., Since residuals of the B1/B2 ionosphere-free code combination) were analyzed as well. Taking the same residuals of the B1/B2 ionosphere-free code combination) were analyzed as well. Taking the same station, for for example, example, Figure Figure 15 15 shows shows the the code code residuals residuals (Res) (Res) of of representative representative IGSO IGSO and and MEO MEO station, satellites on on XMIS. XMIS. As As expected, expected, the the uncorrected uncorrected elevation-dependent elevation-dependent code code bias bias variations variations are are visible visible satellites in the code residuals for both IGSO and MEO satellites. Such variations are stronger for MEOs since in the code residuals for both IGSO and MEO satellites. Such variations are stronger for MEOs since

Sensors 2016, 16, 909

13 of 16

5.3. Code Residual Since the un-modeled code biases will be partly reflected in the residuals, the code residuals (i.e., residuals of the B1/B2 ionosphere-free code combination) were analyzed as well. Taking the same station, for example, Figure 15 shows the code residuals (Res) of representative IGSO and MEO satellites on XMIS. As expected, the uncorrected elevation-dependent code bias variations are visible in the code residuals for both IGSO and MEO satellites. Such variations are stronger for MEOs since they are more seriously affected by the code-phase divergence. After application of the code bias correction models, the residuals run much more like zero mean random noises, which implies that most of the code biases have been well settled. Moreover, residuals of all the involved IGSO and MEO satellites were grouped together and the statistics were shown in Figure 16 and Table 5. The code residuals show a mean of 9 mm for the IGSO group satellites, and a zero mean for the MEO group satellites for all the three schemes. For IGSOs, residual distributions of the latter two schemes (i.e., code biases corrected by the traditional and improved correction models) are slightly better than that of the first scheme (i.e., none code bias corrections applied). As for the MEOs, code residuals show much better distributions in terms of standard deviation after mitigating the code biases with either the traditional or the improved correction models. Sensors 2016, 16, 909 13 of 15

60

0

30

None

0

Traditional

Res Res [m] [m]

-5 5

Improved

0 -5 5 0 -5 8

MEO(C12) Res Res [m] [m]

90

Res Res [m] [m]

Res Res [m] [m]

Res Res [m] [m]

Res Res [m] [m]

C10(IGSO) 5

12

16 UTC [h]

20

24

90

5

60

0

30

None

-5 5

0

Elevation Elevation [deg] [deg]

13 of 15

Elevation Elevation [deg] [deg]

Sensors 2016, 16, 909

0 Traditional

-5 5 0

Improved

-5 16

17

18

19 20 UTC [h]

21

22

Figure Timeseries seriesofofrepresentative representativecode code residuals residuals on Figure 15.15.Time on XMIS XMIS(DoY (DoY130/2015). 130/2015). Traditional

Improved

20

20

15

Percentage Percentage [%] [%]

Percentage Percentage [%] [%]

None 25

15 10 5 0

-5

0 5

-5 0 5 Residual [m]

-5

Traditional

Improved

10 5 0

0 5

None

-5

0 5

-5 0 5 Residual [m]

-5

0 5

Figure Histograms thecode coderesiduals residuals(left: (left: IGSOs; IGSOs; right: Figure 16.16. Histograms ofofallallthe right: MEOs) MEOs)on onXMIS XMIS(DoY (DoY130/2015). 130/2015). Table 5. Mean standard deviation (STD) values of residuals code residuals with different processing Table 5. Mean andand standard deviation (STD) values of code with different processing schemes schemes (unit: m). (unit: m).

Mean STD Mean

None None 0.009 1.166 0.009

STD

1.166

IGSO IGSO Traditional Improved Traditional Improved 0.009 0.009 1.155 1.152 0.009 0.009 1.155

1.152

None None 0.000 1.802 0.000 1.802

MEO MEO Traditional Improved Traditional Improved 0.000 0.000 1.730 1.740 0.000 0.000 1.730

1.740

6. Conclusions Different from the other GNSS systems, BeiDou (IGSO and MEO) code measurements are polluted by the so-called satellite-induced code biases, or code-phase divergences, or code bias variations. The multipath combination (MP), which is a geometry-free and ionosphere-free combination, was used to investigate the code bias variations. It was confirmed that such code bias variations are orbit type-,

Sensors 2016, 16, 909

14 of 16

6. Conclusions Different from the other GNSS systems, BeiDou (IGSO and MEO) code measurements are polluted by the so-called satellite-induced code biases, or code-phase divergences, or code bias variations. The multipath combination (MP), which is a geometry-free and ionosphere-free combination, was used to investigate the code bias variations. It was confirmed that such code bias variations are orbit type-, frequency-, and elevation-dependent. In general, the code measurements of MEO satellites tracked on the B1 frequency at high elevations are more likely to be affected by the code bias variations. To mitigate these unexpected variations, an improved correction model was developed in this paper. Approximately two years data recorded on a globally distributed MGEX stations equipped with different receiver types were used. The continuous piecewise linear functions were employed to produce the correction values as well as their precision indexes. To obtain the best-fitting results, the elevation angles were limited to 5˝ –85˝ with a node separation of 10˝ . To make full use of the improved model, the stochastic information of the corrections were encouraged to take into account. PPP tests were conducted with different code bias variation treatments to verify the effectiveness of the proposed model. Positioning results show that the traditional model does not really yield noticeable improvements. Sometimes it may even degrade the positioning accuracy due to its inaccurate stochastic model. However, the positioning accuracy, as well as convergence, were obviously enhanced by the improved model. This implies that the precision of the uncorrected code measurements and even the corrected code measurements with traditional correction model may be too optimistic. However, a more real and reliable stochastic model can be set up if the precision of corrections is taken into account. Moreover, the wide-lane ambiguities and code residuals were derived for comparison, which further confirmed the effectiveness of the proposed model. Acknowledgments: The authors gratefully acknowledge IGS Multi-GNSS Experiment (MGEX) for providing GNSS data and products. This study would not have been possible without such data. Thanks also go to the National Natural Science Foundation of China (No: 41404006, No: 41474025) and the Open Research Fund of State Key Laboratory of Information Engineering in Survey, Mapping and Remote Sensing, China (No. 15P02). Author Contributions: Fei Guo and Wanke Liu conceived the idea and designed the experiments. Fei Guo and Xin Li performed the experiments and analyzed the data. All authors read, revised and approved the final manuscript. Conflicts of Interest: The authors declare no conflict of interest.

References 1.

2. 3. 4. 5. 6. 7.

8.

Montenbruck, O.; Hauschild, A.; Steigenberger, P.; Hugentobler, U.; Teunissen, P.; Nakamura, S. Initial Assessment of the COMPASS/BeiDou-2 Regional Navigation Satellite System. GPS Solut. 2013, 17, 211–222. [CrossRef] Zhao, Q.; Guo, J.; Li, M.; Qu, L.; Hu, Z.; Shi, C.; Liu, J. Initial results of precise orbit and clock determination for COMPASS navigation satellite system. J. Geod. 2013, 87, 475–486. [CrossRef] Deng, Z.; Zhao, Q.; Springer, T.; Prange, L.; Uhlemann, M. Orbit and Clock Determination-BeiDou. In Proceedings of the International GNSS Service (IGS) Workshop, Pasadena, CA, USA, 23–27 June 2014. Wu, X.; Hu, X.; Wang, G.; Zhong, H.; Tang, C. Evaluation of COMPASS ionospheric model in GNSS positioning. Adv. Space Res. 2013, 51, 959–968. [CrossRef] Zhang, R.; Song, W.; Yao, Y.; Shi, C.; Lou, Y.; Yi, W. Modeling regional ionospheric delay with ground-based BeiDou and GPS observations in China. GPS Solut. 2015, 19, 1–10. [CrossRef] Xu, A.; Xu, Z.; Ge, M.; Xu, X.; Zhu, H.; Sui, X. Estimating zenith tropospheric delays from BeiDou navigation satellite system observations. Sensors 2013, 13, 4514–4526. [CrossRef] [PubMed] Montenbruck, O.; Steigenberger, P.; Hauschild, A. Differential Code Bias Estimation using Multi-GNSS Observations and Global Ionosphere Maps. In Proceedings of the Institute of Navigation (ION) International Technical Meeting, San Diego, CA, USA, 26–28 January 2014. Guo, F.; Zhang, X.; Wang, J. Timing group delay and differential code bias corrections for BeiDou positioning. J. Geod. 2015, 89, 427–445. [CrossRef]

Sensors 2016, 16, 909

9. 10. 11. 12. 13. 14. 15.

16. 17.

18. 19. 20. 21. 22. 23.

24. 25. 26. 27. 28. 29.

30. 31. 32. 33.

15 of 16

Shi, C.; Fan, L.; Li, M.; Liu, Z.; Gu, S.; Zhong, S.; Song, W. An enhanced algorithm to estimate BDS satellite’s differential code biases. J. Geod. 2015, 90, 161–177. [CrossRef] Zhang, B.; Teunissen, P.J.G. Characterization of multi-GNSS between-receiver differential code biases using zero and short baselines. Chin. Sci. Bull. 2015, 21, 1840–1849. [CrossRef] Hauschild, A.; Montenbruck, O.; Sleewaegen, J.M.; Huisman, L.; Teunissen, P. Characterization of compass M-1 signals. GPS Solut. 2012, 16, 117–126. [CrossRef] Wang, G.; de Jong, K.; Zhao, Q.; Hu, Z.; Guo, J. Multipath analysis of code measurements for BeiDou geostationary satellites. GPS Solut. 2015, 19, 129–139. [CrossRef] Tang, W.; Deng, C.; Shi, C.; Liu, J. Triple-frequency carrier ambiguity resolution for Beidou navigation satellite system. GPS Solut. 2014, 18, 335–344. [CrossRef] Zhang, X.; He, X. Performance analysis of triple-frequency ambiguity resolution with BeiDou observations. GPS Solut. 2015, 20, 269–281. [CrossRef] Yang, Y.; Li, J.; Wang, A.; Xu, J.; He, H.; Guo, H.; Dai, X. Preliminary assessment of the navigation and positioning performance of BeiDou regional navigation satellite system. Sci. China Earth Sci. 2014, 57, 144–152. [CrossRef] Li, M.; Qu, L.; Zhao, Q.; Guo, J.; Su, X.; Li, X. Precise Point Positioning with the BeiDou Navigation Satellite System. Sensors 2014, 14, 927–943. [CrossRef] [PubMed] Li, X.; Zhang, X.; Ren, X.; Fritsche, M.; Wickert, J.; Schuh, H. Precise positioning with current multi-constellation Global Navigation Satellite Systems: GPS, GLONASS, Galileo and BeiDou. Sci. Rep. 2015, 5. [CrossRef] [PubMed] Li, X.; Ge, M.; Dai, X.; Ren, X.; Fritsche, M.; Wickert, J.; Schuh, H. Accuracy and reliability of multi-GNSS real-time precise positioning: GPS, GLONASS, BeiDou, and Galileo. J. Geod. 2015, 89, 607–635. [CrossRef] Guo, F.; Zhang, X.; Wang, J.; Ren, X. Modeling and assessment of triple-frequency BDS precise point positioning. J. Geod. 2016. [CrossRef] Nadarajah, N.; Teunissen, P.J.G.; Raziq, N. BeiDou inter-satellite-type bias evaluation and calibration for mixed receiver attitude determination. Sensors 2013, 7, 9435–9463. [CrossRef] [PubMed] Nadarajah, N.; Teunissen, P.J.G.; Sleewaegen, J.M.; Montenbruck, O. The mixed-receiver BeiDou inter-satellite-type bias and its impact on RTK positioning. GPS Solut. 2014, 19, 357–368. [CrossRef] Wanninger, L.; Beer, S. BeiDou satellite-induced code pseudorange variations: Diagnosis and therapy. GPS Solut. 2015, 19, 639–648. [CrossRef] Montenbruck, O.; Hauschild, A.; Steigenberger, P.; Hugentobler, U.; Riley, S. A COMPASS for Asia: First Experience with the BeiDou-2 Regional Navigation System. Poster at IGS Workshop in Olsztyn, Poland 2012. Available online: http://www.igs.org/presents/poland2012/posters (accessed on 10 May 2015). Montenbruck, O.; Rizos, C.; Weber, R.; Weber, G.; Neilan, R.; Hugentobler, U. Getting a grip on multi-GNSS-the international GNSS service MGEX campaign. GPS World 2013, 24, 44–49. De Bakker, P.F.; Tiberius, C.C.J.M.; van der Marel, H.; van Bree, R.J.P. Short and zero baseline analysis of GPS L1 C/A, L5Q, GIOVE E1B, and E5aQ signals. GPS Solut. 2012, 16, 53–64. [CrossRef] Springer, T.; Dilssner, F. SVN49 and Other GPS Anomalies; Inside GNSS: Eugene, OR, USA, 2009; pp. 32–36. Montenbruck, O.; Hauschild, A.; Steigenberger, P.; Langley, R.B. Three’s the challenge. GPS World 2010, 21, 8–19. Pearson, K. Note on regression and inheritance in the case of two parents. Proc. R. Soc. Lond. 1895, 58, 240–242. [CrossRef] Montenbruck, O.; Steigenberger, P.; Khachikyan, R.; Weber, G.; Langley, R.B.; Mervart, L.; Hugentobler, U. IGS-MGEX: Preparing the Ground for Multi-Constellation GNSS Science; Inside GNSS: Eugene, OR, USA, 2014; Volume 9, pp. 42–49. Dow, J.M.; Neilan, R.E.; Rizos, C. The International GNSS Service in a changing landscape of Global Navigation Satellite Systems. J. Geod. 2009, 83, 191–198. [CrossRef] Zhang, X.; Andersen, O.B. Surface ice flow velocity and tide retrieval of the Amery ice shelf using precise point positioning. J. Geod. 2006, 80, 171–176. [CrossRef] Boehm, J.; Heinkelmann, R.; Schuh, H. A global model of pressure and temperature for geodetic applications. J. Geod. 2007, 81, 679–683. [CrossRef] Petit, G.; Luzum, B. IERS Conventions 2010 (IERS Technical Note No. 36); International Earth Rotation and Reference Systems Service: Frankfurt, Germany, 2010.

Sensors 2016, 16, 909

34.

35. 36.

37.

16 of 16

Rizos, C.; Montenbruck, O.; Weber, R.; Neilan, R.; Hugentobler, U. The IGS MGEX Experiment as a Milestone for a Comprehensive Multi-GNSS Service. In Proceedings of Institute of Navigation Pacific Positioning, Navigation and Timing (ION Pacific PNT), Honolulu, HI, USA, 23–25 April 2013; pp. 289–295. Wu, J.; Wu, S.; Hajj, G.; Bertiger, W.; Lichten, S. Effects of antenna orientation on GPS carrier phase. Manuscr. Geod. 1993, 18, 91–98. Melbourne, W.G. The case for ranging in GPS-based geodetic systems. In Proceedings of the first International Symposium on Precise Positioning with the Global Positioning System, Rockville, MD, USA, 15–19 April 1985; pp. 373–386. Wübbena, G. Software developments for geodetic positioning with GPS using TI 4100 code and carrier measurements. In Proceedings of the 1st International Symposium on Precise Positioning with the Global Positioning Systems, Rockville, MD, USA, 15–19 April 1985; pp. 403–412. © 2016 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC-BY) license (http://creativecommons.org/licenses/by/4.0/).