2. IONOSPHERE-FREE COMBINATION MODEL

The ionosphere-free combination model can be expressed as follows (Hoffman-Wellenhof et al., Reference Hofmann-Wellenhof, Lichtenegger and Collins2013; Xu, Reference Xu2007):

(1) $$\eqalign{\lambda _{IF}\nabla \Delta \phi _{IF} & = \displaystyle{{\,f_1^2 } \over {\,f_1^2 - f_2^2 }}\lambda _1\nabla \Delta \phi _1 - \displaystyle{{\,f_2^2 } \over {\,f_1^2 - f_2^2 }}\lambda _2\nabla \Delta \phi _2 \cr & = (\nabla \Delta \rho + \nabla \Delta O + \nabla \Delta T + \nabla \Delta M) + \displaystyle{{\,f_1^2 } \over {\,f_1^2 - f_2^2 }}\lambda _1\nabla \Delta N_1 - \displaystyle{{\,f_2^2 } \over {\,f_1^2 - f_2^2 }}\lambda _2\nabla \Delta N_2} $$

Where $\nabla \Delta $ is the double difference operator between the receiver and satellite; λ _IF is the wavelength of the ionosphere-free observation $\nabla \Delta \varphi _{IF} $ ; $\nabla \Delta \varphi _1 $ and $\nabla \Delta \varphi _2 $ are the double difference L1 and L2 carrier phase observables; $\nabla \Delta \rho $ is the double difference satellite-earth range; $\nabla {\it \Delta} O$ is the double difference orbit error; $\nabla \Delta T$ is the double difference tropospheric delay; $\nabla {\it \Delta} M$ is the double difference multipath error and f ₁, f ₂, λ ₁, λ ₂ and $\nabla \Delta N_1 $ , $\nabla {\it \Delta} N_{\rm 2} $ are the frequencies, wavelengths and ambiguities of carrier phase L1 and L2 respectively.

Considering $\nabla \Delta N_2 = \nabla \Delta N_1 - \nabla \Delta N_W $ , then Equation (1) can be rewritten as (Hoffman-Wellenhof et al., Reference Hofmann-Wellenhof, Lichtenegger and Collins2013; Xu, Reference Xu2007):

(2) $$\lambda _W \nabla \Delta \varphi _{IF} = \nabla \Delta \rho + \nabla \Delta O + \nabla \Delta T + \nabla \Delta M + \displaystyle{{\,f_2} \over {\,f_1 - f_2}} \lambda _N \nabla \Delta N_W - \lambda _N \nabla \Delta N_1 $$

Where λ _W  = 86·2 cm, λ_n = 10·7 cm, are the wavelength of the wide-lane and narrow-lane observation respectively, $\nabla \Delta N_W $ is the double difference wide-lane ambiguity, which can be fixed by the MW(Melbourne-Wubbena) combination. As the satellite-earth range between the satellite and the reference station can be calculated precisely, as well as the orbit error and multipath error, the main influence of the determination of $\nabla \Delta N_1 $ in Equation (2) comes from the double difference residual of tropospheric delay. If the length of the baseline is short, the tropospheric delays at the two ends of the baseline do not appear to be much different and can be calculated by the model directly. Thus, the double difference residual of the tropospheric delay becomes small and can be ignored. However, when the length of the baseline is long, the double difference residual of the tropospheric delay should be considered because the tropospheric delays at the two ends of the baseline vary due to the low spatial correlation. Meanwhile, as the influence of the troposphere is related to the elevation angle of the satellite, a mapping function is introduced. Thus, the tropospheric delay can be expressed as the product of the zenith tropospheric delay Ztd and the mapping function MF(E), which is the function of the elevation angle of the satellite. In fact, the zenith tropospheric delay Ztd is composed of the dry and wet component, which are zenith hydrostatic delay (Zhd) and zenith wet delay (Zwd), respectively. The dry component Zhd is estimated through the Global Pressure and Temperature (GPT) model while the wet component Zwd is estimated as unknown parameters in the observation model. As for the mapping function, the Hopfield model is adopted here.

Hence, $\nabla \Delta T_{AB}^{ij} $ , the double difference tropospheric delay between baseline AB and satellite pair i-j can be described as follows (Hoffman-Wellenhof et al., Reference Hofmann-Wellenhof, Lichtenegger and Collins2013; Xu, Reference Xu2007):

(3) $$\eqalign{\nabla \Delta T_{AB}^{ij} & = \nabla \Delta T_{wet} + \nabla \Delta T_{dry} \cr & = Zwd_A[MF(E_A^i ) - MF(E_A^j )] - Zwd_B[MF(E_B^i ) - MF(E_B^j )] + \nabla \Delta T_{dry}} $$

$\nabla \Delta T_{wet} $ , $\nabla \Delta T_{dry} $ are the double difference wet and dry component of the tropospheric delay. After the fixing of the wide-lane ambiguity, considering Equation (2) as well as ignoring the influence of the satellite orbit error and multipath error, the observation equation can be shown as:

(4) $$\eqalign{&\lambda _W \nabla {\it \Delta} \phi _{IF} - \nabla {\it \Delta} \rho - \displaystyle{{\,f_2} \over {\,f_1 - f_2}} \lambda _N \nabla {\it \Delta} N_W - \nabla {\it \Delta} T_{dry} \cr & = \left[ {\matrix{ {\nabla MF(E_A^{ij} )} & { - \nabla MF(E_B^{ij} )} & {\lambda _N} \cr}} \right]\left[ {\matrix{ {Zwd_A} \cr {Zwd_B} \cr {\nabla {\it \Delta} N_1} \cr}} \right]}$$

$\nabla $ is the single-difference operator between the satellite, then the error equation can be written:

(5) $$V = \left[ {\matrix{ {\nabla MF(E_A^{ij} )} & { - \nabla MF(E_B^{ij} )} & {\lambda _N} \cr}} \right]\left[ {\matrix{ {Zwd_A} \cr {Zwd_B} \cr {\nabla \Delta N_1} \cr}} \right] - L$$

where

(6) $$L = \lambda _W \nabla \Delta \phi _{IF} - \nabla \Delta \rho - \displaystyle{{\,f_2} \over {\,f_1 - f_2}} \lambda _N \nabla \Delta N_W - \nabla \Delta T_{dry} $$

Therefore, Equation (5) can be solved by the method of the sequential least square adjustment or Kalman filter algorithm. In the process of the sequential least square adjustment, the observation data from one minute or longer is suggested for the first adjustment, as the observation amount of the first epoch is less than the estimated parameters. In addition, the correlation of the observation data from the neighbouring epochs is quite strong when the interval of the collected data is low. As for the Kalman filter algorithm, the basic ionosphere-free observation equation is shown in Equation (4). Thus, when the two reference stations have n common view satellites (not including the reference satellite), the Kalman filter model can be established as:

(7) $$\left\{ \matrix{\mathop {X_k} \limits_{n + 2,1} = \mathop {X_{k - 1}} \limits_{n + 2,1} + \mathop {\Omega _k} \limits_{n + 2,1} \matrix{ {} & {} & {} & {} & {\Omega _k \sim N(0,D_\Omega )} \cr} \hfill \cr \mathop {L_k} \limits_{n,1} = \mathop {B_k} \limits_{n,n + 2} \cdot \mathop {X_k} \limits_{n + 2,1} + \mathop {\Delta _n} \limits_{n,1} \matrix{ {} & {} & {} & {} & {\Delta _n \sim N(0,D_\Delta )} \cr} \hfill} \right.$$

Where X is the state vector, as $X = [\matrix{ {Zwd_A} & {Zwd_B} & {\nabla {\it \Delta} N_1} & \cdots & {\nabla {\it \Delta} N_n} \cr} ]^T $ ; D _Ω is the variance-covariance matrix of the dynamic noise Ω _k . Due to the slow change of the zenith tropospheric delay in a short time, its noise variance $\sigma _T^2 $ is set from 0·0001 to 0·0009 m²/h as the random walking parameter while the noise variance of the ambiguity $\sigma _N^2 $ is set as zero. D _Δ is the variance-covariance matrix of the observation which can be set by the weight of the carrier phase observation. Further, the weight can be decided through the elevation angle method and so on (Gao et al., Reference Gao, Zhao and Wang2005; He and Yang, Reference He and Yang2001). Thus, the zenith tropospheric delay, float ambiguity solution as well as the corresponding covariance matrix can be estimated in real-time during the filter process. Then, the integer ambiguity can be fixed using the Least-squares Ambiguity Decorrelation Adjustment (LAMBDA) method (Teunissen, Reference Teunissen1995).

The ordinary model of the ionosphere-free combination mentioned above can achieve excellent performance for satellites with high elevation angles. But the success rate of ambiguity resolution for the satellites with low elevation angle is low due to the observation noise and ill-posed problems of the equation.

3. NETWORK RTK AMBIGUITY FIXING MODEL WITH A PRIORI ATMOSPHERIC INFORMATION RESTRICTED

To improve the success rate of the ambiguity resolution, the ambiguities of the satellites can be divided into two categories; one kind is the ambiguities that can be fixed easily as N _e , the other is the opposite ones that are difficult to fix as N _h , thus Equation (5) can be written as:

(8) $$\left( {\matrix{ {V_1} \cr {V_2} \cr}} \right) = \left( {\matrix{ {MF_1} & {C_1} & O \cr {MF_2} & O & {C_2} \cr}} \right)\left( {\matrix{ {T_{zwd}} \cr {N_e} \cr {N_h} \cr}} \right) - \left( {\matrix{ {L_1} \cr {L_2} \cr}} \right)$$

Where V ₁ and V ₂ are the residuals of the error equations for solving the ambiguity N _e and N _h , respectively. C ₁ and C ₂ are the wavelengths of the narrow-lane combination, respectively. $T_{zwd} = [\matrix{ {Zwd_A} & {Zwd_B} } ]^T $ , $MF_{} = [\matrix{ {\nabla MF(E_A^{ij} )} & { - \nabla MF(E_B^{ij} )}} ]$ . Based on the least squares principle, if the two categories are calculated together, the consequence of the estimated parameters is as follows:

(9) $$\eqalign{&\left( {\matrix{ {\hat T_{zwd}} \cr {\hat N_e} \cr {\hat N_h} \cr}} \right) = \left( {\matrix{ {MF_1^T P_1 MF_1 + MF_2^T P_2 MF_2} & {MF_1^T P_1 C_1} & {MF_2^T P_2 C_2} \cr {C_1^T P_1 MF_1} & {C_1^T P_1 C_1} & O \cr {C_2^T P_2 MF_2} & O & {C_2^T P_2 C_2} \cr}} \right)^{ - 1} \cr &\hskip4.3pc\quad\times\left( {\matrix{ {MF_1^T P_1 L_1 + MF_2^T P_2 L_2} \cr {C_1^T P_1 L_1} \cr {C_2^T P_2 L_2} \cr}} \right)}$$

In the process above, the advantage of the satellite ambiguities that can be fixed easily is not taken. In view of this, the satellite ambiguities that can be fixed easily can be fixed first, and then the fixed ambiguities can be brought into the equation as known parameters. Thus Equation (8) can be transformed as:

(10) $$\left( {\matrix{ {V_1} \cr {V_2} \cr}} \right) = \left( {\matrix{ {MF_1} & O \cr {MF_2} & {C_2} \cr}} \right)\left( {\matrix{ {T_{zwd}} \cr {N_h} \cr}} \right) - \left( {\matrix{ {L_1 - C_1 N_e} \cr {L_2} \cr}} \right)$$

The consequence of the estimation is:

(11) $$\eqalign{&\left( {\matrix{ {\hat T_{zwd}} \cr {\hat N_h} \cr}} \right) = \left( {\matrix{ {MF_1^T P_1 MF_1 + MF_2^T P_2 MF_2} & {MF_1^T P_1 C_2} \cr {C_2^T P_2 MF_2} & {C_2^T P_2 C_2} \cr}} \right)^{ - 1} \cr &\hskip4.5pc\times \left( {\matrix{ {MF_1^T P_1 (L_1 - C_1 N_e ) + MF_2^T P_2 L_2} \cr {C_2^T P_2 L_2} \cr}} \right)}$$

In this model, the structure of the model becomes stronger with fewer parameters as well as the T _zwd restriction condition. Furthermore, the model can achieve more accurate float ambiguity resolution as the ill-posed problem reduces.

4. DESIGN OF THE AMBIGUITY FIXING SCHEME IN THE INITIALISATION OF THE NETWORK RTK SYSTEM

In the process of the initialisation of the Network RTK, different satellite ambiguities with different elevation angles show different difficulties in dealing with the ambiguity resolution. For satellites with high elevation angles, their atmospheric delay is small and can be corrected completely by the atmospheric model. Thus, their ambiguities can be fixed using several epochs or just one epoch without the estimation of the zenith tropospheric delay. In this case, if the double difference atmospheric delay obtained from the easily fixed satellites can be used to restrict the unfixed satellite ambiguity resolution, a rapid fixing solution can be achieved. In addition, to distinguish the high and low elevation angle, 35° is selected as the boundary in this paper. Though there is no strict standard to define which elevation angle is high or low, the boundary of 30° to 35° is often accepted (Deng and Wang, Reference Deng and Wang2015). It is well known that the lower the elevation angle, the larger the mapping function. Hence, a small bias of the zenith tropospheric delay will result in a large difference of the tropospheric delay along the line of sight, which will affect the resolution of the ambiguities. In view of this, 35° is chosen in this paper. Therefore, the three-step rapid ambiguity fixing algorithm with a priori atmospheric information restricted can be adopted for the first initialisation of the Network RTK system. The specific calculation scheme is shown in Figure 1.

Figure 1. Three-step ambiguity fixing scheme for the reference stations within the Network RTK.

5. DATA ANALYSIS.

5.1. Case 1: Data from Jiangsu Continuously Operating Reference Stations (JSCORS)

To validate the actual effect of the proposed method, real data collected from two stations (NJTS and NJPK) of the JSCORS network from China, which form a 51 km baseline, are employed. Data were collected on 29 December 2012, from 0:00:00 to 1:00:00, with a collection of 1000 epochs. The observed satellites are: G01, G04, G17, G20, G28, G32 for GPS and R01, R08, R14, R15, R17, R24 for GLONASS. The reference satellites are G28 and R24, respectively. The elevations of the satellites at station NJTS are shown in Figure 2.

Figure 2. Elevation angle at station NJTS.

Based on the fixed wide-lane ambiguity $\nabla \Delta N_1 $ , the double difference L1 carrier phase ambiguity can be correctly fixed with the use of the ionosphere-free combination model. In this case, the satellite pairs with high elevation angles are: G01-G28, G17-G28, G20-G28, G32-G28, R01-R24 and R14-R24, for their elevation angles are above 35°.

For these satellites, the influence of the troposphere can be ignored as the delay becomes small after being amended by the tropospheric model as well as double differencing. Therefore, their ambiguities can be solved within one single epoch as ambiguities become the sole unknown parameters in Equation (4). As can be seen from Figure 3, the double difference ambiguities of high elevation angle satellites can actually be fixed well within one single epoch. Furthermore, with regard to satellites with low elevation angles, the ambiguity resolution performance of the proposed algorithm is presented below, which is compared with the conventional method.

Figure 3. Ambiguity difference between float and fixed solution within single epoch for high elevation angle satellites.

In this case, the sequential least square adjustment is selected for the three-step algorithm. Figure 4(a) shows the performance of the conventional method. As we can see, for some satellites, the difference between float solution and fixed solution are kept at about two cycles after an initialisation of about 200 epochs. Combined with the ratio value shown in Figure 5, it is evident that the ratio value remains at a low level, meaning the reliability of the ambiguity resolution is poor. The performance of the algorithm proposed in this paper is shown in Figure 4(b). It is observed that, for the ambiguity resolution of the satellites with low elevation angles, the ambiguity difference between float and fixed solution is less than one cycle after the first epoch. Besides, Figure 5 shows that the ratio value rises up at first, indicating a high level reliability. Meanwhile, from Figure 4(b) we can see that there exist some fluctuations for the float solution. The main reason is that the obtained tropospheric delay, which is used to restrict the ambiguity resolution of the low elevation angle satellite, has some differences with the actual tropospheric delay at the coordinate position of the low satellites due to the influence of the mapping function. In addition, from comparison of the condition number of the normal equation of the two methods, as shown in Figure 6, we can conclude that the condition number of the normal equation of the conventional method reduces as the observation time grows. However, when applied to the proposed algorithm, the condition number of the normal equation is several orders of magnitude smaller than the conventional method, and varies little. It indicates that the proposed algorithm can reduce the ill-posed problems of the model, and improve the speed of the ambiguity resolution.

Figure 4. Comparison of the difference of the float and fixed ambiguity solution using the two methods.

Figure 5. Ratio value of the two methods.

Figure 6. Condition number of the normal equation.

5.2. Case 2: Data from Suzhou Continuously Operating Reference Stations (SZCORS)

To validate the general effect of the proposed algorithm further, data from three reference stations, jint, taoy and weit, are selected from SZCORS, part of the JSCORS network. The distribution of the stations is shown in Figure 7. The lengths of the baselines, jint-taoy, taoy-weit and weit-jin, are 40 km, 52 km and 66 km. In this case, an arbitrary continuous 2000 epochs are selected and the interval is one second. In this period, there were 15 GPS/GLONASS satellites in common view. The distribution of the satellites in the sky is shown in Figure 8.

Figure 7. Relative distribution of the stations.

Figure 8. Satellite distribution in the sky.

According to the difference of the elevation angle, the satellites are divided into two categories; one is the high elevation angle satellites and the other is the low elevation angle satellites, as shown in Table 1.

Table 1. The two categories of the satellites.

Based on the correctly fixed wide-lane ambiguities, ionosphere-free observations are constructed. The conventional method and the proposed algorithm are both applied to calculate the GPS/GLONASS double difference ambiguities of the three baselines. Then we can calculate the coordinates of the stations by substituting the ambiguities calculated every epoch into Equation (12). In fact, the coordinates of the reference stations are precisely known. Therefore we can acquire the bias of the calculated coordinates based on the precisely known coordinates of the reference stations. In this way, we can check the correctness of the calculated ambiguities through the coordinate bias.

(12) $$\eqalign{& \lambda _w \nabla {\it \Delta} \phi _{IF} - \displaystyle{{\,f_1} \over {\,f_1 + f_2}} \nabla {\it \Delta} N_1 \cdot \lambda _w + \displaystyle{{\,f_2} \over {\,f_1 + f_2}} \nabla {\it \Delta} N_2 \cdot \lambda _w - (\nabla \rho _A^{\,jr} - \nabla \rho _{B0}^{\,jr} ) - \nabla \Delta T_{dry} \cr & = [\matrix{ {\nabla l_B^{\,jr}} & {\nabla m_B^{\,jr}} & {\nabla n_B^{\,jr}} & {\nabla MF(E_A^{ij} )} & { - \nabla MF(E_B^{ij} )} \cr} ] \cdot \left( {\matrix{ {\delta x} \cr {\delta y} \cr {\delta z} \cr {Zwd_A} \cr {Zwd_B} \cr}} \right)}$$

Where $\nabla l_B^{\,jr} $ , $\nabla m_B^{\,jr} $ and $\nabla n_B^{\,jr} $ are the coefficients of the linearized double difference satellite-earth range of coordinate component x, y and z. The coordinate bias of the North (N) and East (E) of the three baselines, jint-taoy, taoy-weit and weit-jint are shown in Figures 9, 10 and 11, where (a) is the result of the conventional method and (b) is the result from the proposed algorithm.

Figure 9. Baseline:jint-taoy.

Figure 10. Baseline: taoy-weit.

Figure 11. Baseline: weit – jint.

The coordinate bias reflects the effect of the methods for ambiguity resolution. From Figure 9 to Figure 11, taking the convergence speed and reliability of the coordinate into consideration, we can see that the performance of the proposed method is better than the conventional one, especially in the convergence speed. In the conventional method, the convergence time is several hundreds or thousands of epochs while the proposed method shows a much lower convergence time. That is to say, the a priori atmospheric information restricted algorithm can obtain the fixed ambiguity quickly. Moreover, comparing the condition number of the normal equation of the two methods which are shown in Figure 12, we can conclude that the proposed algorithm can reduce the ill-posed problems of the NRTK model effectively and thus improve the speed of the ambiguity resolution.

Figure 12. Condition number of the normal equation.