1. INTRODUCTION
The BeiDou Navigation Satellite System (BDS) has been officially providing Positioning, Navigation and Timing (PNT) services for the Asia-Pacific region since 27 December 2012 (Shi et al., Reference Shi, Zhao, Hu and Liu2013). By the end of 2014, the system comprised 14 satellites available for services, all broadcasting triple-frequency signals. The centre frequency values are 1,561·098 MHz (B1), 1,207·140 MHz (B2) and 1,268·520 MHz (B3) (Yang et al., Reference Yang, Li, Xu, Tang, Guo and He2011; Li et al., Reference Li, Zhang, Ren, Fritsche, Wickert and Schuh2015a). Multi-frequency signals not only contribute more code and carrier phase observables, but also bring many benefits in terms of the efficiency and reliability of Ambiguity Resolution (AR), which are crucial for real-time precise positioning applications (Li et al., Reference Li, Ge, Zhang, Nischan and Wickert2013a; Reference Li, Ge, Dai, Ren, Fritsche, Wickert and Schuh2015b). Since the European Global Navigation Satellite System-2 (GNSS-2) (later called Galileo) and American modernised Global Positioning System (GPS) were proposed in the last century, there have been many significant research efforts into the carrier phase AR using triple-frequency signals. The earlier studies, such as the Three-Carrier Ambiguity Resolution (TCAR) method (Forssell et al., Reference Forssell, Martin-Neira and Harrisz1997; Vollath et al., Reference Vollath, Birnbach and Landau1998) and the Cascade Integer Resolution (CIR) method (De Jonge et al., Reference De Jonge, Teunissen, Jonkman and Joosten2000; Hatch et al., Reference Hatch, Jung, Enge and Pervan2000), essentially used the same geometry-free model to estimate the carrier phase ambiguity with three-step rounding procedures (Teunissen et al., Reference Teunissen, Joosten and Tiberius2002). It is a requirement for reliable AR that the total noise level of the virtual code and carrier phase must be distinctly smaller than the wavelength. This requirement can be easily satisfied for Extra-Wide-Lane (EWL, λ ≥ 2·93 m) and Wide-Lane (WL, 2·93 > λ ≥ 0·75 m) combinations, but it is hardly complied with for the Medium-Lane (ML, 0·75 > λ ≥ 0·19 m) and Narrow-Lane (NL, λ < 0·19 m) combinations, particularly in the long-range case (Li et al, Reference Li, Feng and Shen2010).
The classical TCAR method employs a geometry-free model to eliminate the geometric-related terms, and then resolves the ambiguity by averaging and rounding operations. The reliability of AR is biased mainly by the ionospheric delay (Jung et al., Reference Jung, Enge and Pervan2000). Feng (Reference Feng2008) presented a general modelling strategy for AR to minimise the effect of the ionospheric delay, which introduced ionosphere-reduced virtual signals to eliminate the ionosphere parameters in the geometry-based observational model for the long baseline case. The strategy outlined three major steps and utilised a geometry-based estimator at steps 2 and 3 to improve the AR success probability. The TCAR concept and algorithm were extended by Feng and Li (Reference Feng and Li2008; Reference Feng and Li2009) in recent studies for introducing both the geometry-based and geometry-free model for AR. Li et al. (Reference Li, Feng and Shen2010) identified that tropospheric delay effects were the key limitation on the geometry-based model TCAR for long baselines, and employed a Geometry-free and Ionosphere-Free (GIF) model to avoid the effect of tropospheric delay on ML/NL observable AR in the long-range case.
Since there were insufficient triple-frequency GPS or Galileo satellites available for earlier studies, performances were mainly validated by theoretical analysis or using a semi-simulation method proposed by Li (Reference Li2008) to generate the third frequency signal based on the existing dual-frequency GPS measurements. For BDS, Tang, et al. (Reference Tang, Deng, Shi and Liu2014) assumed that the standard deviation of the B3 code observation noise is significantly smaller than that of B1 and B2, and introduced a geometry-based model stepwise AR algorithm, which determined float ambiguities using the standard least-squares adjustment and resolved the carrier phase integer ambiguities using the Least squares Ambiguity Decorrelation Adjustment (LAMBDA) method for each step. Zhang and He (Reference Zhang and He2015) used real BeiDou observation data to test the performances of LAMBDA, geometry-free TCAR, geometry-based TCAR and geometry-free and ionospheric-free TCAR methods, and concluded that the geometry-based TCAR and LAMBDA methods have good AR performance for short baselines. However, the tropospheric delay effect, which cannot be eliminated by the geometry-based TCAR and LAMBDA methods, is still the major problem for long-range cases. Wang and Rothacher (Reference Wang and Rothacher2013) derived a geometry-free and ionosphere-free combination for GPS, Galileo and COMPASSIII and analysed the characteristics of this combined method. Zhao et al. (Reference Zhao, Dai, Hu, Sun, Shi and Liu2015) employed a modified TCAR method to eliminate or reduce the residual ionospheric delay in steps 2 and 3 of the classical TCAR, but it did not take the different characteristics of triple-frequency code observation noises into consideration, and did not show an encouraging improvement in WL and NL success rates.
In this paper, we mainly focus on the geometry-free TCAR method and introduce an improved TCAR method that is suitable for BDS, which can overcome the ionospheric and tropospheric effects to provide reliable AR. The characteristics of the BeiDou triple-frequency code observation noises have been analysed and the limitations on the geometry-free TCAR method for providing high reliability and probability of AR have been taken into consideration, particularly for the long baseline case. The performance of the improved TCAR method has been evaluated by real BeiDou observation data for one short baseline and two long baseline cases. The research findings of this paper will be summarised in the final section.
2. GEOMETRY-FREE MODEL FOR TCAR
2.1. Fundamental observation equations and definitions
The centre frequencies of the B1, B2 and B3 bands are set as f 1, f 2 and f 3, respectively. Taking the impacts of the satellite orbital error, the tropospheric delay and the first-order ionospheric delay into consideration, the fundamental Double-Differenced (DD) code and carrier phase observation equations can be expressed as:
$$\Delta P_i = \Delta \rho + \Delta \delta _{orb} + \Delta \delta _{tro} + \displaystyle{{\,f_1^2} \over {\,f_i^2}} \Delta \delta I_1 + \varepsilon _{\Delta Pi} $$
$$\Delta \Phi _i = \Delta \rho + \Delta \delta _{orb} + \Delta \delta _{tro} - \displaystyle{{\,f_1^2} \over {\,f_i^2}} \Delta \delta I_1 - \lambda _i \cdot \Delta N_i + \varepsilon _{\Delta \Phi i} $$
where the symbol “Δ” represents the DD operator product; ΔP i and ∆Φ i are the DD code and carrier phase measurements for frequency f i in metres; Δρ is the DD geometric distance between the satellite and the receiver; Δδ orb denotes the DD orbital error; Δδ tro is the DD tropospheric delay; ΔδI 1 represents the DD first-order ionospheric delay for frequency f 1; and ε∆Pi and ε∆Φi are the DD code and carrier phase observation noises for frequency f i , respectively.
The combined triple-frequency DD carrier phase observation can be expressed as (Feng, Reference Feng2008):
$$\matrix{ {\Delta \Phi _{\left( {i,j,k} \right)}} \hfill & { = \displaystyle{{i \cdot f_1 \cdot \Delta \Phi _1 + j \cdot f_2 \cdot \Delta \Phi _2 + k \cdot f_3 \cdot \Delta \Phi _3} \over {i \cdot f_1 + j \cdot f_2 + k \cdot f_3}}} \hfill \cr \hfill & { = \Delta \rho + \Delta \delta _{orb} + \Delta \delta _{tro} - \beta _{\left( {i,j,k} \right)} \cdot \Delta \delta I_1 - \lambda _{\left( {i,j,k} \right)} \cdot \Delta N_{\left( {i,j,k} \right)} + \varepsilon _{\Delta \Phi \left( {i,j,k} \right)}} \hfill \cr} $$
where the combined coefficients i, j, k are all arbitrary integers. The combined wavelength is defined as:
$$\lambda _{\left( {i,j,k} \right)} = \displaystyle{c \over {i \cdot f_1 + j \cdot f_2 + k \cdot f_3}} = \displaystyle{{\lambda _1 \cdot \lambda _2 \cdot \lambda _3} \over {i \cdot \lambda _2 \cdot \lambda _3 + j \cdot \lambda _1 \cdot \lambda _3 + k \cdot \lambda _1 \cdot \lambda _2}} $$
The combined integer ambiguity is:
$$\Delta N_{\left( {i,j,k} \right)} = i \cdot \Delta N_1 + j \cdot \Delta N_2 + k \cdot \Delta N_3 $$
The combined first-order ionospheric scale factor is defined as:
$$\beta _{\left( {i,j,k} \right)} = \displaystyle{{\,f_1^2 \left( {{i / {\,f_1 + {\,j / {\,f_2}} + {k / {\,f_3}}}}} \right)} \over {i \cdot f_1 + j \cdot f_2 + k \cdot f_3}} $$
In contrast with the definition of the carrier phase combination, we denote the combined triple-frequency DD code observation as:
$$\eqalign{{\Delta {P_{\left( {l,m,n} \right)}} \vert l + m + n = 1} & = {l \cdot \Delta {P_1} + m \cdot \Delta {P_2} + n \cdot \Delta {P_3}} \cr & = {\Delta \rho + \Delta {\delta _{orb}} + \Delta {\delta _{tro}} + {\beta _{\left( {l,m,n} \right)}} \cdot \Delta \delta {I_1} + {\varepsilon _{\Delta P\left( {l,m,n} \right)}}}} $$
The combined first-order ionospheric scale factor for the code observation is:
$$\beta _{\left( {l,m,n} \right)} = l + \displaystyle{{\,f_1^2} \over {\,f_2^2}} m + \displaystyle{{\,f_1^2} \over {\,f_3^2}} n = q_1 \cdot l + q_2 \cdot m + q_3 \cdot n$$
where q i = f 1 2/f i 2 is the ratio of the ionospheric delay effects on frequency f i to that on frequency f 1.
2.2. General formation of the geometry-free TCAR model
The term “geometry-free model” refers to an observational model for the ambiguity parameter that is formed without any geometric-related term (Feng and Rizos, Reference Feng and Rizos2009), such as the geometric distance from the satellite to the receiver, the orbital error or the tropospheric delay, etc. For this type of model, the carrier phase integer ambiguity is generally estimated by rounding the float value to its nearest integer. A general geometry-free observational model for a virtual integer parameter can be expressed as (Feng et al., Reference Feng, Rizos and Higgins2007):
$$\Delta N_{\left( {i,j,k} \right)} = \displaystyle{{\Delta P_{\left( {l,m,n} \right)} - \Delta \Phi _{\left( {i,j,k} \right)}} \over {\lambda _{\left( {i,j,k} \right)}}} - \displaystyle{{\beta _{\left( {l,m,n} \right)} + \beta _{\left( {i,j,k} \right)}} \over {\lambda _{\left( {i,j,k} \right)}}} \Delta \delta I_1 - \displaystyle{{\varepsilon _{\Delta P\left( {l,m,n} \right)} - \varepsilon _{\Delta \Phi \left( {i,j,k} \right)}} \over {\lambda _{\left( {i,j,k} \right)}}} $$
The accuracy of the geometry-free model AR will be affected by two types of errors: the first is the type of random errors, which includes the combined carrier phase and the code observation noise, and generally obeys a normal distribution with a zero mean parameter. The second is the type of systematic errors, which is generally affected by the combined ionospheric delay. This type of error has a non-zero mean value and cannot be eliminated or reduced by averaging over a long period of epochs. The effects of these two types of errors for AR reliability are different and should be distinguished. A weighted noise level in cycles is introduced to identify the effects of different errors for AR purposes in the geometry-free model:
$$\sigma _{WTC} = \displaystyle{{\sqrt {\gamma \cdot \left( {\beta _{\left( {l,m,n} \right)} + \beta _{\left( {i,j,k} \right)}} \right)^2 \cdot \left( {\Delta \delta I_1} \right)^2 + \varepsilon _{\Delta P\left( {l,m,n} \right)}^2 + \varepsilon _{\Delta \Phi \left( {i,j,k} \right)}^2}} \over {\lambda _{\left( {i,j,k} \right)}}} $$
where γ is the weighted factor for the first-order ionospheric delay. The value of γ is correlated with the values of the wavelength, ionospheric delay and observation noise, which is extraordinarily complicated and will be discussed in the following section of this paper. Generally, the weighted factor meets the condition γ ≥1.
2.3. The real time estimation method of the BDS triple-frequency code observation noise
Since the carrier phase noise is fairly small, usually considered to be 0·01 cycles for practical purposes, we can assume that the carrier phase noises for each of the three frequencies are independent and have an identical standard deviation, i.e., ε ∆Φ1 = ε ∆Φ2 = ε ∆Φ3 ≡ε ∆Φ. The linearly combined carrier phase observation noise can be expressed as (Li et al., Reference Li, Feng and Shen2010; Zhao et al., Reference Zhao, Dai, Hu, Sun, Shi and Liu2015):
$$\varepsilon _{\Delta \Phi \left( {i,j,k} \right)}^2 = \displaystyle{{\left( {i \cdot f_1} \right)^2 + \left( {\,j \cdot f_2} \right)^2 + \left( {k \cdot f_3} \right)^2} \over {\left( {i \cdot f_1 + j \cdot f_2 + k \cdot f_3} \right)^2}} \cdot \varepsilon _{\Delta \Phi} ^2 \equiv \mu _{\left( {i,j,k} \right)}^2 \cdot \varepsilon _{\Delta \Phi} ^2 $$
Since the code chipping rates, lock loop parameters, and multipath effects are distinctly different for each of the BeiDou triple-frequency code observables, the code observation noise level for each frequency is diverse and requires a real-time estimation. The noise levels and multipath effects for the carrier phase observable are much smaller than those for the code observable and can be ignored (Wang et al., Reference Wang, de Jong, Zhao, Hu and Guo2015). Therefore, the code observation noise can be analysed using a combination of code and carrier phase observables by eliminating the geometric-related terms and the first-order ionospheric delay effects:
$$\varepsilon _{{\Delta P} _i} = \Delta P_i - \left( {c_{1i} \cdot \Delta \Phi _1 + c_{2i} \cdot \Delta \Phi _2 + c_{3i} \cdot \Delta \Phi _3} \right)$$
where c ki (k = 1, 2, 3) represents the combined coefficient. The mathematical expression of the constraint condition is:
$${\rm M}in\left\{ {\left. {c_{1i}^2 + c_{2i}^2 + c_{3i}^2} \right \vert c_{1i} + c_{2i} + c_{3i} = 1;\displaystyle{1 \over {\,f_i^2}} + \displaystyle{{c_{1i}} \over {\,f_1^2}} + \displaystyle{{c_{2i}} \over {\,f_2^2}} + \displaystyle{{c_{3i}} \over {\,f_3^2}} = 0} \right\}$$
By utilising the BDS triple-frequency values to solve the constraint functions, the code observation noise series can be approximately extracted as follows:
$$\eqalign{& \varepsilon _{{\Delta P} _1} = \Delta P_1 - \left( {4 \cdot 1665 \cdot \Delta \Phi _1 - 2 \cdot 3483 \cdot \Delta \Phi _2 - 0 \cdot 8182 \cdot \Delta \Phi _3} \right) \cr & \varepsilon _{{\Delta P} _2} = \Delta P_2 - \left( {5 \cdot 2424 \cdot \Delta \Phi _1 - 3 \cdot 1010 \cdot \Delta \Phi _2 - 1 \cdot 1414 \cdot \Delta \Phi _3} \right) \cr & \varepsilon _{{\Delta P} _3} = \Delta P_3 - \left( {4 \cdot 9897 \cdot \Delta \Phi _1 - 2 \cdot 9242 \cdot \Delta \Phi _2 - 1 \cdot 0655 \cdot \Delta \Phi _3} \right)} $$
The extracted noise series has constant bias terms, such as the unresolved carrier phase integer ambiguities and the receiver hardware delays. Assuming that there are no cycle slips over observing periods, the bias terms can be eliminated through averaging over a long period of epochs. The standard deviation of DD code observation noise can be estimated as follows in real time:
$$\varepsilon _{\Delta Pi} \left( K \right) = \sqrt {\displaystyle{1 \over {K - 1}}\sum\limits_{n = 1}^K {\left( {\varepsilon _{\Delta Pi} \left( n \right) - \displaystyle{{\sum\limits_{n = 1}^K {\varepsilon _{\Delta Pi} \left( n \right)}} \over K}} \right)^2}} $$
In order to illustrate the noise levels on DD code observables for each frequency, we collected BDS three-day triple-frequency data from the Day of Year (DOY) 280 to 282, 2014 at two stations with 179 m distances in Changsha, China. The satellite C03, which can be continuously visible and at higher elevation angle, was set as pivot satellite, and the elevation cut-off angle was set to 15°. The standard deviations of the DD code observation noise are shown in Table 1.
Table 1. Standard deviations of the BDS DD code observation noise.
Since the code chipping rates, the modulation recognition techniques of the signals and the lock loop parameters of the receiver are different for each frequency, the standard deviations of the DD code observation noises are frequency-dependent and satellite-dependent. As can be observed in Table 1, since the chipping rate on B3 is ten times faster than that on B1 and B2, the standard deviation of ∆P 3 is smaller than that for ∆P 1 and ∆P 2. However this is inconsistent with the assumption proposed by Tang et al. (Reference Tang, Deng, Shi and Liu2014) that the standard deviation of ∆P 3 is practically 0·2 (theoretically it is 0·1) times than that of ∆P 1 and ∆P 2. So in practical application, it is reasonable and necessary to estimate the noise levels of the triple-frequency code observables in real time.
3. OPTIMAL COMBINATION OF CARRIER PHASE AND CODE OBSERVABLES
The traditional TCAR method seeks to determine the suitable combinations of carrier phase and code measurements by searching and comparing operations within the integer range. Therefore this method cannot always achieve a minimal total noise level, and also requires a large number of computations and comparisons, limiting the AR success probability.
3.1. Optimal combined coefficients of code observables
Given a group of definite combined coefficients i, j, k for the carrier phase observables, and since the combined wavelength λ (i, j, k) is positive, the constraint condition to minimise the weighted noise level σ W TC can be expressed as follows:
$$Min\left\{ {\left. {F\left( {l,m,n} \right) = \gamma \cdot \left( {\beta _{\left( {l,m,n} \right)} + \beta _{\left( {i,j,k} \right)}} \right)^2 \cdot \left( {\Delta \delta I_1} \right)^2 + \varepsilon _{\Delta P\left( {l,m,n} \right)}^2 + \varepsilon _{\Delta \Phi \left( {i,j,k} \right)}^2 \vert l + m + n = 1} \right\}} \right.$$
This problem is therefore equivalent to finding the minimum value in the constraint function. We introduce the Lagrange function f (l, m, n, θ) as follows:
$$\eqalign{{\,f\left( {l,m,n,\theta } \right)} &= {\gamma \cdot {{\left( {{\beta _{\left( {l,m,n} \right)}} + {\beta _{\left( {i,j,k} \right)}}} \right)}^2} \cdot {{\left( {\Delta \delta {I_1}} \right)}^2} + {l^2} \cdot \varepsilon _{\Delta P1}^2 + {m^2} \cdot \varepsilon _{\Delta P2}^2 + {n^2} \cdot \varepsilon _{\Delta P3}^2 } \cr & \quad + {\mu _{\left( {i,j,k} \right)}^2 \cdot \varepsilon _{\Delta \Phi }^2 + \theta \left( {l + m + n - 1} \right)}} $$
Calculating the first-order partial derivative of the independent variables and setting them to zero gives:
$$\eqalign{& \displaystyle{{\partial f\left( {l,m,n,\theta} \right)} \over {\partial l}} = 2q_1 \cdot \gamma \cdot \left( {q_1 \cdot l + q_2 \cdot m + q_3 \cdot n + \beta _{\left( {i,j,k} \right)}} \right) \cdot \left( {\Delta \delta I_1} \right)^2 +\; 2l \cdot \varepsilon _{\Delta P1}^2 + \theta = 0 \cr & \displaystyle{{\partial f\left( {l,m,n,\theta} \right)} \over {\partial m}} = 2q_2 \cdot \gamma \cdot \left( {q_1 \cdot l + q_2 \cdot m + q_3 \cdot n + \beta _{\left( {i,j,k} \right)}} \right) \cdot \left( {\Delta \delta I_1} \right)^2 +\; 2m \cdot \varepsilon _{\Delta P2}^2 + \theta = 0 \cr & \displaystyle{{\partial f\left( {l,m,n,\theta} \right)} \over {\partial n}} = 2q_3 \cdot \gamma \cdot \left( {q_1 \cdot l + q_2 \cdot m + q_3 \cdot n + \beta _{\left( {i,j,k} \right)}} \right) \cdot \left( {\Delta \delta I_1} \right)^2 +\; 2n \cdot \varepsilon _{\Delta P3}^2 + \theta = 0 \cr & \displaystyle{{\partial f\left( {l,m,n,\theta} \right)} \over {\partial \theta}} = l + m + n - 1 = 0} $$
Solving the equations, the optimal combined coefficients l, m, n for the code observables can be calculated through the following equations:
$$\eqalign{&l = \displaystyle{{\matrix{\varepsilon _{\Delta P2}^2 \cdot \varepsilon _{\Delta P3}^2 + \gamma \left( {q_3^2 - {q_1}{q_3} + \left( {{q_3} - {q_1}} \right) \cdot {\beta _{\left( {i,j,k} \right)}}} \right){{\left( {\Delta \delta {I_1}} \right)}^2} \cdot \varepsilon _{\Delta P2}^2 \cr +\; \gamma \left( {q_2^2 - {q_1}{q_2} + \left( {{q_2} - {q_1}} \right) \cdot {\beta _{\left( {i,j,k} \right)}}} \right){{\left( {\Delta \delta {I_1}} \right)}^2} \cdot \varepsilon _{\Delta P3}^2 }} \over {\matrix{\varepsilon _{\Delta P1}^2 \cdot \varepsilon _{\Delta P2}^2 + \varepsilon _{\Delta P1}^2 \cdot \varepsilon _{\Delta P3}^2 + \varepsilon _{\Delta P2}^2 \cdot \varepsilon _{\Delta P3}^2 \cr +\; \gamma \cdot \left( {{{\left( {{q_2} - {q_3}} \right)}^2}\varepsilon _{\Delta P1}^2 + {{\left( {{q_1} - {q_3}} \right)}^2}\varepsilon _{\Delta P2}^2 + {{\left( {{q_1} - {q_2}} \right)}^2}\varepsilon _{\Delta P3}^2 } \right) \cdot {{\left( {\Delta \delta {I_1}} \right)}^2}}}} \cr & {m = \displaystyle{{\matrix{\varepsilon _{\Delta P1}^2 \cdot \varepsilon _{\Delta P3}^2 + \gamma \left( {q_3^2 - {q_2}{q_3} + \left( {{q_3} - {q_2}} \right) \cdot {\beta _{\left( {i,j,k} \right)}}} \right)\cr {{\left( {\Delta \delta {I_1}} \right)}^2} \cdot \varepsilon _{\Delta P1}^2 + \gamma \left( {q_1^2 - {q_1}{q_2} + \left( {{q_1} - {q_2}} \right) \cdot {\beta _{\left( {i,j,k} \right)}}} \right){{\left( {\Delta \delta {I_1}} \right)}^2} \cdot \varepsilon _{\Delta P3}^2 }} \over {\matrix{\varepsilon _{\Delta P1}^2 \cdot \varepsilon _{\Delta P2}^2 + \varepsilon _{\Delta P1}^2 \cdot \varepsilon _{\Delta P3}^2 + \varepsilon _{\Delta P2}^2 \cdot \varepsilon _{\Delta P3}^2 \cr + \;\gamma \cdot \left( {{{\left( {{q_2} - {q_3}} \right)}^2}\varepsilon _{\Delta P1}^2 + {{\left( {{q_1} - {q_3}} \right)}^2}\varepsilon _{\Delta P2}^2 + {{\left( {{q_1} - {q_2}} \right)}^2}\varepsilon _{\Delta P3}^2 } \right) \cdot {{\left( {\Delta \delta {I_1}} \right)}^2}}}}} \cr & {n = \displaystyle{{\matrix{\varepsilon _{\Delta P1}^2 \cdot \varepsilon _{\Delta P2}^2 + \gamma \left( {q_2^2 - {q_2}{q_3} + \left( {{q_2} - {q_3}} \right) \cdot {\beta _{\left( {i,j,k} \right)}}} \right)\cr{{\left( {\Delta \delta {I_1}} \right)}^2} \cdot \varepsilon _{\Delta P1}^2 + \gamma \left( {q_1^2 - {q_1}{q_3} + \left( {{q_1} - {q_3}} \right) \cdot {\beta _{\left( {i,j,k} \right)}}} \right) {{\left( {\Delta \delta {I_1}} \right)}^2} \cdot \varepsilon _{\Delta P2}^2 }} \over {\matrix{\varepsilon _{\Delta P1}^2 \cdot \varepsilon _{\Delta P2}^2 + \varepsilon _{\Delta P1}^2 \cdot \varepsilon _{\Delta P3}^2 + \varepsilon _{\Delta P2}^2 \cdot \varepsilon _{\Delta P3}^2 \cr + \;\gamma \cdot \left( {{{\left( {{q_2} - {q_3}} \right)}^2}\varepsilon _{\Delta P1}^2 + {{\left( {{q_1} - {q_3}} \right)}^2}\varepsilon _{\Delta P2}^2 + {{\left( {{q_1} - {q_2}} \right)}^2}\varepsilon _{\Delta P3}^2 } \right) \cdot {{\left( {\Delta \delta {I_1}} \right)}^2}}}}}}$$
For BDS, the ratios q 1, q 2 and q 3 have values of 1, 1·67 and 1·51, respectively. Therefore, given a group of combined coefficients for the carrier phase observables, the optimal combined coefficients for the code observables are mainly affected by the following factors: the code observation noise for each frequency, the weighted factor and the first-order ionospheric delay. Since we have already introduced a real time estimation method for the BDS triple-frequency code observation noise, the other two factors will be discussed in the following sections.
3.2. Estimation of the ionospheric delay
The ionospheric delay effect is a major error source that limits the AR success rates of the geometry-free model. In particular, for the long baseline case, as well as the reduction on spatial correlation, the DD ionospheric delay significantly increases, by up to several metres, and seriously jeopardises precise positioning applications (Li et al., Reference Li, Zhang and Ge2011; Reference Li, Ge, Zhang and Wickert2013b).
Assuming that the two EWL/WL observables ΔΦ(i1,j1,k1) and ΔΦ(i2,j2,k2) have been obtained and their integer ambiguities ΔN (i1,j1,k1) and ΔN (i2,j2,k2) have been resolved, a simple method of estimation of DD ionospheric delay can be expressed as follows (Feng and Rizos, Reference Feng and Rizos2005):
$$\Delta \delta \hat I_1 = \displaystyle{{\matrix{\Delta \Phi _{\left( {i1,j1,k1} \right)} - \Delta \Phi _{\left( {i2,j2,k2} \right)} - \lambda _{\left( {i1,j1,k1} \right)} \cdot \Delta N_{\left( {i1,j1,k1} \right)} \cr + \lambda _{\left( {i2,j2,k2} \right)} \cdot \Delta N_{\left( {i2,j2,k2} \right)} - \varepsilon _{\Delta \Phi \left( {i1,j1,k1} \right)} + \varepsilon _{\Delta \Phi \left( {i2,j2,k2} \right)} }} \over {{\beta _{\left( {i2,j2,k2} \right)} - \beta _{\left( {i1,j1,k1} \right)}}}} $$
However, the observation noises are amplified by the combination operation, this estimation method for ionospheric delay is extremely noisy and the standard deviation of this estimator is over 60εΔΦ. A smoothing operation can improve the estimated accuracy, but will lose the specific variation between consecutive epochs. To overcome this problem, Hatch filtering is employed:
$$\Delta \delta {\breve I_1}\left( t \right) = \left( {1 - \omega } \right) \cdot \left({\Delta \delta {{\breve I} _1}} \left( {t - 1} \right) + \Delta \delta {{\tilde I}_1} {\left( {t,t - 1} \right)}\right) + \omega \cdot \Delta \delta {{\breve I} _1} \left( t \right)$$
where Δδ
${{\breve I} _1} $
(t) is the estimated value of the first order ionospheric delay after Hatch filtering at time t; and the symbol ω is the smoothing factor and can be set as the 1/n, where n is the epoch number that requires smoothing. Δδ
${{\tilde I} _1} $
(t,t-1) is the variation of the ionospheric delay over the observation interval, and can be accurately calculated by obtaining the difference between the B1 and B3 DD carrier phase observables over two consecutive epochs to eliminate the geometric-related terms and carrier phase integer ambiguities as follows:
$$\Delta \delta \tilde I_1 \left( {t,t - 1} \right) = \displaystyle{ \eqalign{&\Delta \Phi _1 \left( t \right) - \Delta \Phi _1 \left( {t - 1} \right) - \left({ \Delta \Phi _3 \left( t \right) - \Delta \Phi _3 \left( {t - 1} \right) }\right)\cr & + \Delta \varepsilon _{\Phi 1} \left( {t,t - 1} \right) - \Delta \varepsilon _{\Phi 3} \left( {t,t - 1} \right)} \over {q_3 - 1}}$$
The standard deviation of the estimation for Δδ
${{\tilde I} _1} $
(t,t-1) is 3·9εΔΦ. So the standard deviation of the Hatch filtering ionospheric delay estimator will be approximately less than 6εΔΦ, which indicates that the estimation accuracy of the ionospheric delay has been significantly improved.
3.3. Discussion on the preferences of the weighted factor
The estimation methods for the code observation noise and ionospheric delay have been introduced in the previous sections. In this section, we will discuss the preferences of the weighted factor for ionospheric delay.
As discussed in the previous analysis, the accuracy of the geometry-free model AR is affected by two types of errors. Both types of errors can be normalised based on the wavelength. The random error is presumed to obey a normal distribution N 1 (0, σ 1 2). Its bias term is zero, and the variance σ 1 2 is equal to the variance of the combined code and carrier phase observation noises. The systematic error can also be considered to have a normal distribution N 2 (μ 2, σ 2 2). Its bias term μ 2 is determined by the estimation of ionospheric delay, and the variance σ 2 2 is equal to the estimated variance. Since the weighted total noise level is the cumulative effect of these two types of errors, it also obeys a normal distribution.
To obtain the geometry-free model integer ambiguity solutions from Equation (9), the direct method is to round the float ambiguity value to its nearest integer value. The success rate can be computed through the integral of normal probability density function (Feng, Reference Feng2008):
$$P\left( { - 0 \cdot 5 \lt x \lt 0 \cdot 5} \right) = \mathop{\int}\limits_{ - 0 \cdot 5}^{0 \cdot 5} {\displaystyle{1 \over {\sigma \sqrt {2\pi}}}} \exp \left( { - \displaystyle{{\left( {x - bias} \right)^2} \over {2\sigma ^2}}} \right)dx$$
Figure 1 illustrates the estimated probability of AR when the float ambiguity value obeys a normal distribution. The bias term is mainly determined by the estimated value of ionospheric delay. Since the variance of combined code observation noise is much larger than that of the combined carrier phase observation noise and the estimated variance for the ionospheric delay, the variance σ 2 is mostly determined by the combined code observation noise.
Figure 1. Estimated probability of AR success rate for normal distribution.
It is obvious that when the absolute value of the bias term is smaller than 0·2 cycles, the effect of the bias term for the estimated probability is approximately equivalent to the variance. However, when it is larger than 0·2 cycles, and in particular larger than 0·5 cycles, the effect of the bias term on the estimated probability significantly increases. In this paper, we present a simple method for the preferences of the weighted factor:
$$\left\{ {\matrix{ {\gamma = 1} \hfill &, \hfill & {\matrix{ {\left \vert {bias} \right \vert \le 0 \cdot 2} & {cycles} \cr}} \hfill \cr {\gamma = round\left(\displaystyle{{\left \vert {bias} \right \vert} \over {0 \cdot 2}}\right)} \hfill &, \hfill & {\matrix{ {\left \vert {bias} \right \vert \gt 0 \cdot 2} & {cycles} \cr}} \hfill \cr}} \right.$$
It should be noted that this method is only a representative example, and there are many other methods that can be used to obtain the preferences of the weighted factor, which can be investigated in future research.
Since the methods of estimation and preferences for the parameters in Equation (19) have been discussed, Table 2 presents pairs of EWL/WL observables and their weighted total noise levels for different baseline cases. Based on the approximate average values in Table 1, we assume that the standard deviations of the BDS triple-frequency code observation noises are: ε ∆ P1 = 0·45 m, ε ∆ P2 = 0·60 m and ε ∆ P3 = 0·30 m; the carrier phase measurement noises for BDS three frequencies are identical and are set as 0·005 m; the DD ionospheric delays are 0·1, 0·4 and 1 m for short, medium and long distances, respectively (Feng, Reference Feng2008); and the weighted factor γ is determined by Equation (21). The signal ∆P 3, the high precision code measurement on B3 (Tang et al., Reference Tang, Deng, Shi and Liu2014), is involved for comparison.
Table 2. The weighted noise level for several EWL/WL observables.
It is obvious that the weighted total noise level is significantly reduced by using the optimal combination. This can be particularly seen for the medium and long range cases, where the weighted total noise level was reduced by approximately 30% for most combinations shown in Table 2. It can also be seen that the first EWL signal ΔΦ(0,−1,1) has the lowest noise level for all baseline cases, which indicates that this observable can be regarded as the optimal combination signal for AR. The other EWL ΔΦ(−1,−5,6) signal has the longest wavelength and good AR performance for the short baseline case, but performs poorly as the length of baseline increases. The WL ΔΦ(1,1,−2) signal has a lower weighted total noise level in long baseline case, and can be selected as the other combination signal for AR in long distances cases.
3.4. Geometry-free and ionosphere-eliminated method for reliable AR of ML/NL observables
In order to compute the basic ambiguities, a minimum of three linearly independent combined observables are required. Therefore, the third combination should be chosen from a new category, which is linearly independent of the previous two EWL/WL observables. The problem is that there are no such EWL/WL candidates in the new category that will allow the ambiguities to be fixed as easily as the previous two observables using the virtual code measurements. Feng and Rizos (Reference Feng and Rizos2005) simply chose the carrier phase measurement ΔΦ(1,0,0) as the third candidate, and smoothed the float values over multiple epochs to fix the ambiguity. Li et al. (Reference Li, Feng and Shen2010) formed a GIF combination as the new category candidate as follows:
$$\eqalign{\Delta N_{\left( {i,j,k} \right)} & = \displaystyle{{\alpha _1 \Delta \tilde \Phi _{\left( {i1,j1,k1} \right)} + \alpha _2 \Delta \tilde \Phi _{\left( {i2,j2,k2} \right)} - \Delta \Phi _{\left( {i,j,k} \right)}} \over {\lambda _{\left( {i,j,k} \right)}}} \cr & \quad - \displaystyle{{\alpha _1 \varepsilon _{\Delta \tilde \Phi \left( {i1,j1,k1} \right)} + \alpha _2 \varepsilon _{\Delta \tilde \Phi \left( {i2,j2,k2} \right)} - \varepsilon _{\Delta \Phi \left( {i,j,k} \right)}} \over {\lambda _{\left( {i,j,k} \right)}}}}$$
where
$\Delta \tilde \Phi _{\left( {i1,j1,k1} \right)} $
and
$\Delta \tilde \Phi _{\left( {i2,j2,k2} \right)} $
are two EWL/WL observables, and their ambiguities have been fixed. Keeping within the constraints of being both geometry-free and ionosphere-free, the combined coefficients a
1 and a
2 are determined by:
$$\left. {\matrix{ {\alpha _1 \beta _{\left( {i1,j1,k1} \right)} + \alpha _2 \beta _{\left( {i2,j2,k2} \right)} - \beta _{\left( {i,j,k} \right)} = 0} \cr {\alpha _1 + \alpha _2 = 1} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\cr}} \right\}$$
The standard deviation in cycles of the ML/NL observable obtained using the GIF method can be expressed as:
$$\sigma \left( {\Delta N_{\left( {i,j,k} \right)}} \right) = \displaystyle{{\sqrt {\alpha _1^2 \varepsilon _{\Delta \tilde \Phi \left( {i1,j1,k1} \right)}^2 + \alpha _2^2 \varepsilon _{\Delta \tilde \Phi \left( {i2,j2,k2} \right)}^2 + \varepsilon _{\Delta \Phi \left( {i,j,k} \right)}^2}} \over {\lambda _{\left( {i,j,k} \right)}}} $$
Although the GIF method is free of the ionospheric and tropospheric effects and can be used without a constraint on distance, it is very noisy for AR on single epoch measurements. Since we have introduced an accurate unbiased estimator for the ionospheric delay, the new category candidates can be found by eliminating the ionospheric delay with the estimated value Δδ
${{\breve I} _1} $
. We investigate an optimal combination which is Geometry-Free and Ionosphere-Eliminated (GFIE) for the new category candidate AR:
$$\eqalign{\Delta {N_{\left( {i,j,k} \right)}} & = {\displaystyle{{{\alpha _1}\Delta {{\tilde \Phi }_{\left( {i1,j1,k1} \right)}} + {\alpha _2}\Delta {{\tilde \Phi }_{\left( {i2,j2,k2} \right)}} - \Delta {\Phi _{\left( {i,j,k} \right)}}} \over {{\lambda _{\left( {i,j,k} \right)}}}}} \cr & + \displaystyle{{{\alpha _1}{\beta _{\left( {i1,j1,k1} \right)}} + {\alpha _2}{\beta _{\left( {i2,j2,k2} \right)}} - {\beta _{\left( {i,j,k} \right)}}} \over {{\lambda _{\left( {i,j,k} \right)}}}} {\Delta \delta } {{\breve I}_1}{\left( t \right)} \cr & - \displaystyle{{{\alpha _1}{\varepsilon _{\Delta \tilde \Phi \left( {i1,j1,k1} \right)}} + {\alpha _2}{\varepsilon _{\Delta \tilde \Phi \left( {i2,j2,k2} \right)}} - {\varepsilon _{\Delta \Phi \left( {i,j,k} \right)}}} \over {{\lambda _{\left( {i,j,k} \right)}}}}}$$
The constraint conditions for GFIE are such that the combination is free of geometric-related terms and the combined noise is minimised:
$$\left. \matrix{ \eqalign{Min \Big\{ \alpha _1^2 \varepsilon_{\Delta \tilde \Phi \left( {i1,j1,k1} \right)}^2 + \alpha _2^2 \varepsilon _{\Delta \tilde \Phi \left({i2,j2,k2}\right)}^2 + \varepsilon _{\Delta \Phi \left({i,j,k}\right)}^2 + (\alpha _1 \beta _{(i1,j1,k1)} \cr + {\alpha _2} \beta _{\left( {i2,j2,k2} \right)} - \beta _{\left({i,j,k} \right)})^2 {\sigma ^2} \left(\Delta \delta {\breve I}_1\right) \Big\} } \cr \hfill {a_1} + {a_2} = 1 \cr } \right\}$$
The Lagrange function can be set up to find the optimal coefficients for the GFIE combination:
$$\matrix{{\,f\left( {{\alpha _1},{\alpha _2},\theta } \right)} = {\alpha _1^2 \varepsilon _{\Delta \tilde \Phi \left( {i1,j1,k1} \right)}^2 + \alpha _2^2 \varepsilon _{\Delta \tilde \Phi \left( {i2,j2,k2} \right)}^2 + \varepsilon _{\Delta \Phi \left( {i,j,k} \right)}^2 } \hfill \cr \quad\quad\quad\quad\quad\quad{ + {{({\alpha _1}{\beta _{\left( {i1,j1,k1} \right)}} + {\alpha _2}{\beta _{\left( {i2,j2,k2} \right)}} - {\beta {\left( {i,j,k} \right)}})}^2}{\sigma ^2}\left( {\Delta \delta ({{\breve I}_1})} \right) + \theta \left( {{a_1} + {a_2} - 1} \right)} \cr } $$
From the previous analysis, it is known that the standard deviation of the estimation for ionospheric delay is less than 6εΔΦ, Table 3 shows the standard deviations in cycles of the GIF and GFIE methods for reliable AR of the new category of ML/NL observables. The carrier phase noise is assumed to be 0·005 m, and the two basic EWL/WL observables are
$\Delta \tilde \Phi _{\left( {0, - 1,1} \right)} $
and
$\Delta \tilde \Phi _{\left( {1,1, - 2} \right)} $
.
Table 3. Performances of GIF and GFIE methods for reliable AR of the new category of ML/NL observables.
It is obvious that the GFIE method not only eliminates the effects of the ionospheric delay and tropospheric delay, but also improves the estimated accuracy for the reliable AR of the new category of ML/NL observables. This makes it possible to achieve a high single epoch AR success rate in geometry-free model without the distance constraint.
4. EXPERIMENT AND ANALYSIS
In this section, we evaluate the performance of the improved TCAR method for BDS. A short baseline experiment with 179 m distances was carried out in Changsha, China from DOY 280 to 281, 2014. Observation data were collected at 30 s sampling interval by two BeiDou triple-frequency receivers manufactured by the Satellite Navigation R&D Center. The choke-ring BDS antennas were located on the roofs of buildings to reduce the multipath effects. The long baseline performance was analysed with the observation data collected from DOY 143 to 144, 2015 at Xi'an, Kunming and Wuhan stations. Both code and carrier phase observation data were available for three frequencies with 30 s sampling interval, and the receivers were manufactured by the 20th Institute of China Electronics Technology Group Corporation (CETC). The Xi'an station was taken as the master station, and the baseline distances were 654 km and 1,205 km to the rover stations of Wuhan and Kunming, respectively.
Figure 2 shows the ambiguity residuals of combination signals ∆N (0,−1,1), ∆N (−1,−5,6) and ∆N (1,1,−2) for the satellite pair C01-C03 in the short baseline case with 179 m distance (left), and for the satellite pair C07-C03 in the long baseline case with 1,205 km distance (right). Compared with the classical code observable ∆P 3, the optimal combination decreases the standard deviation of float ambiguities, especially for the long range case. Since the optimal combination reduces the ionospheric delay effects, the float ambiguities are closer to the integer value with less bias components, which is more conducive to the geometry-free model AR. It can also be found that the EWL signal ∆N (−1,−5,6) shows good performance for the short baseline case, however the float ambiguities can be observed to follow the ionospheric delay changes for the long-range case. The WL signal ∆N (1,1,−2) shows good performance for both short and long baseline conditions, especially for long distance; the bias component is mitigated by using the optimal combination. Although the standard deviations are still large to achieve perfect (100%) AR success probability for single epoch measurements, they benefit from the distribution of float ambiguities, and a perfect success probability can be efficiently achieved by averaging the AR results over several epochs in the time domain.
Figure 2. Geometry-free model TCAR performances of the optimal code combination (blue point) and the classical combination ∆P 3 (red plus sign) for the AR of the EWL observables ∆N (0,−1,1) (top), ∆N (−1,−5,6) (middle) and the WL observable ∆N (1,1,−2) (bottom) for the satellite pair C01-C03 for the short baseline case with 179 m distance (left), and for the satellite pair C07-C03 for the long baseline case with 1,205 km distance (right).
A statistical analysis algorithm is used to evaluate the improved TCAR performance, and the AR success probability can be defined as the percentage of the correctly solved epoch number to the total epoch number:
$$P\left( s \right) = \displaystyle{{N_{suc}} \over {N_{all}}} \times 100\% $$
where P(s) is the AR success probability, N suc is the epoch number for successful resolution and N all is the total number of epochs. The settings for the statistical algorithm are listed as follows: the elevation angles of visible satellites should be higher than 15° and the C03 is set as the pivot satellite. Table 4 shows the success probability of several EWL/WL observables in different range cases for the BDS service. The classical code measurement ∆P 3 is involved for comparison. The results show that both the optimal and classical code combinations display a good performance for the AR of EWL/WL observables for a short baseline. When the length of baseline increases, the AR success probability significantly decreases for the classical combination, whereas in contrast, the optimal combination still displays good performance. The success probabilities increase by at least 10% for most EWL/WL signals by using the optimal combination.
Table 4. Success probability of geometry-free model AR for several EWL/WL observables.
The efficiencies of the GIF and GFIE AR methods will now be assessed for a new category of ML/NL candidates. As illustrated in Figure 3, the GFIE method has remarkably reduced the standard deviation of the float ambiguities for the ML/NL observables, in both short and long baseline cases. Thus, it becomes efficient to fix the integer ambiguities by using the rounding operation.
Figure 3. The performance of GIF method (red plus sign) and GFIE method (blue point) for ∆N (2,0,−3) (top) and ∆N (0,0,1) (bottom) AR for the satellite pair C01-C03 in the short baseline case with 179 m distance (left), and for the satellite pair C07-C03 in the long baseline case with 1,205 km distance (right).
Table 5 lists the single epoch AR success probability for the new category of ML/NL observables. Compared with the GIF method, the GFIE method significantly improves the single epoch AR success probability from below 24% to above 80%. Since the GFIE method is also free of both ionospheric delay and geometric-related effects, the success probability can be further improved by obtaining an average over several epochs.
Table 5. Single epoch AR success probability of GIF and GFIE methods for the new category of ML/NL observables.
5. CONLUSION AND OUTLOOK
In this paper, we have analysed the limitations of the geometry-free TCAR method. Since the noise levels of BDS triple-frequency code observables are distinctly different, and the ionospheric delay majorly impacts the geometry-free model AR, we have introduced estimation methods for the noise level of code observable and DD ionospheric delay in real time, and have proposed an improved strategy for the geometry-free TCAR method which is suitable for BDS service. This strategy introduces a method for finding the optimal combination of code observables to minimise the total noise level, and also investigates the GFIE method to obtain high AR success probability for the third method linearly independent of ML/NL observables.
To verify and demonstrate the proposed strategy, both short and long baseline BeiDou triple-frequency observation data have been collected and processed to analyse the performance of the improve geometry-free TCAR method. Compared with the traditional geometry-free TCAR method, the improved TCAR method exhibits stable and reliable performance for all baselines. The optimal combination reduces the total noise level and increases the AR success probability by more than 10% for most EWL/WL combinations in the long-range case. The AR success probability for the new category of ML/NL observables is also significantly improved by using the GFIE method, and the single epoch AR success probability is even above 80% for the long baseline case. This result is very promising for long baseline real-time kinematic applications.
