1. INTRODUCTION
With the modernisation of GNSS and the rapid development of satellite systems, more satellite and frequency data can be used for precise positioning. The geometric strength of the observation satellite can be enhanced by the multi-system data combination, and it is beneficial for improving positioning accuracy, reliability and reducing initialisation time (Force and Miller, Reference Force and Miller2013; Odolinski et al., Reference Odolinski, Teunissen and Odijk2014; Tegedor et al., Reference Tegedor, Øvstedal and Vigen2014). Multi-system satellite data has great advantages over single-system data, especially in the complex environment, and it is obviously reflected in relative positioning (Pan et al., Reference Pan, Meng, Gao, Wang and Dodson2014; Teunissen et al., Reference Teunissen, Odolinski and Odijk2014; Liu et al., Reference Liu, Shu, Xu, Qian, Zhang and Zhang2017).
There are mainly two models that were used for relative positioning. The first one is the classic combination in which each system uses its reference satellite to form the double-differenced (DD) model, and then combine the multi-system double-difference observation solution. Which also known as the loose combination model. The other is that all systems select one reference satellite and form the DD model with it. Not only can the double-difference observation equation be formed in a system, but also form the DD equation between systems by the way. This combination is also called the tight combination model (Julien et al., Reference Julien, Alves, Cannon and Zhang2003). The tight combination model can add additional observation equations to the loose combined model, and the positioning is more usable in complex environments (Odijk et al., Reference Odijk, Nadarajah, Zaminpardaz and Teunissen2016).
The influence of DISB should be considered when applying tight combination model. The code DISB can be estimated together with the coordinate, but the ambiguity will absorb the phase DISB and affect the positioning result (Odijk and Teunissen, Reference Odijk and Teunissen2013a). In this case, although the ambiguity of the DD observations in the system still keeps as an integer, whereas the ambiguity of the double-difference observations between systems will not. In fact, only the fractional part of the phase DISB parameter affects the fixation of the ambiguity, because the integer part of the phase DISB can be absorbed by the ambiguity without affecting its integer characteristics. If the fractional part of the phase DISB parameter can be accurately obtained in advance, it will not affect the ambiguity fixed.
For the tight combination model in overlapping frequency, such as the GPS L1 and the Galileo E1, the phase DISB can be calculated by selecting the reference satellites in different satellite systems for parameter recombination, and the fractional part of the result is selected as the phase DISB results finally (Odijk and Teunissen, Reference Odijk and Teunissen2013a, Reference Odijk and Teunissen2013b). For the tight combination model of different frequencies, the phase DISB can also be calculated by the parameter recombination method. However, the calculation results include the deviation of the fractional part caused by the ambiguity transformation, the impact should be considered (Gao et al. Reference Gao, Gao and Pan2017, Reference Gao, Meng and Gao2018). However, the parameter recombination method requires long-term observations to be accurately calculated. In addition, Paziewski and Wielgosz (Reference Paziewski and Wielgosz2015) proposed an approach that uses a priori value to calculate the phase DISB parameter (Paziewski and Wielgosz, Reference Paziewski and Wielgosz2015), but the accuracy of the method depends on the accuracy of the prior value. Tian et al. (Reference Tian, Ge and Neitzel2015) proposed to use the relationship between the phase DISB and the ratio value, use the particle filter to calculate the phase DISB (Tian et al., Reference Tian, Ge and Neitzel2015, Reference Tian, Ge and Neitzel2017), but the method is computationally intensive. In summary, a simple and effective phase DISB calculation method has a certain significance.
From the above analysis, a direct calculation method of phase DISB parameters is proposed based on the correspondence between phase DISB and ratio values in this paper. According to the ratio calculation formula, the method performs a spatial transformation on the ambiguity floating value and directly calculates the phase DISB by using the geometric relationship. Different from the particle filter, the result of the method just relies on a value that can make the ratio maximum. Because a large number of particle calculations are unnecessary, the computation amount is reduced obviously.
The rest of this paper is organised as follows: The second part introduces the calculation formula of DISB parameters, the third part introduces the calculation process of this method in detail, and finally the data verification and summary discussion of the method.
2. DISB PARAMETER CALCULATION METHOD
2.1. Tight combination positioning model
The un-differenced (UD) model of GNSS pseudorange and phase phase observations is expressed as (Teunissen and Kleusberg, Reference Teunissen, Kleusberg, Kleusberg and Teunissen1996):
$$\begin{cases} \varphi _{j,k}^s =\rho _k^s +c(dt_k -dt^s)-I_{j,k}^s +T_k^s +\lambda _j (N_{j,k}^s +\delta _{j,k} +\delta _j^S )+b_k^s +\varepsilon _{j,k}^s \\ P_{j,k}^s =\rho _k^s +c(dt_k -dt^s)+I_{j,k}^s +T_k^s +d_{j,k} +d_j^S +b_k^s+e_{j,k}^s \end{cases}$$
where ϕ is phase observation and P is the pseudorange observation. s, j and k represent the satellite s, frequency j and receiver k. ρ is the distance between satellite and receiver, c is the light speed in vacuum, dt k and dt s are the clock offsets for receiver and satellite. 
$I_{j, k}^{s}$ is the ionospheric delay, 
$T_{k}^{s}$ is the tropospheric delay. λj is the wavelength under the frequency j, 
$N_{j, k}^{s}$ is the ambiguity, δj, k and 
$\delta _{j}^{S}$ are the hardware delay for the receiver and constellation S on phase observation, d j, k and 
$d_{j}^{S}$ are the hardware delay for the receiver and constellation S on pseudorange observation, 
$\varepsilon _{j, k}^{s} $ and 
$e_{j, k}^{s} $ are the measurement noise, respectively.
During the calculation, the satellite clock is calculated by using the broadcast ephemeris, and the errors such as the sagnac and the relativistic effect were corrected by the model. Short baselines are mainly considered in this paper, the atmospheric error and multipath error between stations can be ignored.
The data in GPS L1/Galileo E1/ BDS B1 were selected to form the tight combination model as the example in this paper. The satellite with the highest altitude in the GPS was selected as the reference satellite and named as 1g. The remaining satellites were numbered as m. The form of the DD observation equation in the GPS is:
$$\left\{ {\begin{array}{@{}l} \varphi_{kl}^{1_g m} =\rho_{kl}^{1_g m} +\lambda_{L1} N_{kl}^{1_g m} +\varepsilon_{kl}^{1_g m} \\ P_{kl}^{1_g m} =\rho_{kl}^{1_g m} +e_{kl}^{1_g m} \\ \end{array}} \right.$$It is assumed that satellites of other constellation S can be observed beside the GPS, the frequency of data is j, and the satellite number is r. The system observation data is used to establish a DD observation equation with the reference satellite 1 g. The equation is:
$$\left\{ {\begin{array}{@{}l} \varphi_{kl}^{1_g r} =\rho_{kl}^{1_g r} +\lambda_j (N_{kl}^r +\delta_{kl}^S )-\lambda_{{\rm L}_1 } (N_{kl}^{1_g } +\delta_{kl}^G )+\varepsilon_{kl}^{1_g r} \\ P_{kl}^{1_g r} =\rho_{kl}^{1_g r} +d_{kl}^{GS} +e_{kl}^{1_g r} \\ \end{array}} \right.$$In the Equation 3, the single-difference ambiguity and the phase hardware delay are correlated which makes the ambiguity and hardware delay parameter results cannot be calculated directly. In response to this problem, the current method is to reorganise the parameters. Assume that satellite 1r in the S system is selected as the reference satellite in the system S. Recompose the above observation equations into the following form:
$$\left\{ {\begin{array}{@{}l} \varphi_{kl}^{1_g 1_S } =\rho_{kl}^{1_g 1_S } +\lambda_j (N_{kl}^{1_S } +\delta_{kl}^S )-\lambda_{{\rm L}_1 } (N_{kl}^{1_g } +\delta_{kl}^G )+\varepsilon_{kl}^{1_g 1_S } \\ P_{kl}^{1_g 1_S } =\rho_{kl}^{1_g 1_S } +d_{kl}^{GS} +e_{kl}^{1_g 1_S } \\ \varphi_{kl}^{1_g r} =\rho_{kl}^{1_g r} +\lambda_j N_{kl}^{r1_S } +\lambda_j (N_{kl}^{1_S } +\delta_{kl}^S )-\lambda_{{\rm L}_1 } (N_{kl}^{1_g } +\delta_{kl}^G )+\varepsilon_{kl}^{1_g r} \\ P_{kl}^{1_g r} =\rho_{kl}^{1_g r} +d_{kl}^{GS} +e_{kl}^{1_g r} \end{array}} \right.$$1S is the reference satellite in system S and remain satellites in it are numbered as r. Other parameters are the same as the Equation 1. The phase DISB between GPS and system S is:
$$\lambda _j \delta _{kl}^{GS} =\lambda _j (N_{kl}^{1_r } +\delta _{kl}^S )-\lambda _{{\rm L}_1 } (N_{kl}^{1_g } +\delta _{kl}^G )$$The unknown parameter can be solved under the combination of the DD equation between two constellations. The unknown parameter matrix is:
$$X=[x,y,z,d_{kl}^{G{S}} ,\delta _{kl}^{GS} ,N_{kl}^{1_g m_1 } ,\ldots ,N_{kl}^{1_g m_{n-1} } ,N_{kl}^{1_S r_1 } ,\ldots ,N_{kl}^{1_S r_{m-1} } ]$$ It should be noted that the phase DISB result contains the reference satellite ambiguity, so there is a deviation from the actual value. When the frequencies j and L1 coincide, the deviation is only the phase hardware delay fractional part. When the frequencies j and L1 are different, the deviation includes hardware delay fractional part and fractional part of 
$( \lambda_{{\rm L}_{1} } N_{kl}^{1_{g} } ) /\lambda_{{\rm j}} $. If the reference satellite of the GPS changes, the deviation value will change accordingly, and the latter parameter can be stabilised by means of DD ambiguity transfer (Gao et al., Reference Gao, Gao and Pan2017). Therefore, the fractional part is taken as the phase DISB parameter in general.
2.2. Estimation method of DISB based on the principle of maximum ratio
When the number of observations in the system S is small, the results of the above calculation methods may be unstable. In addition, the result obtained by only one method is not reliable enough. To develop another method of the DISB calculating is meaningful.
In the tight combination model described above, the result of un-reorganised model is:
$$X=[x,y,z,d_{kl}^{GS} ,N_{kl}^{1_g s_1 } ,\ldots ,N_{kl}^{1_g s_{n-1} } ,\bar {\delta }_{kl}^{GS_1} ,\ldots ,\bar {\delta }_{kl}^{G{S}_q } ]$$ Where 
$\bar {\delta }_{kl}^{GS_{r} } =( N_{kl}^{r} +\delta _{kl}^{S} ) -( \lambda _{{\rm L}_{1} } /\lambda _{j} ) ( N_{kl}^{1_{g} } +\delta _{kl}^{G} ) $, n is the number of observations in the GPS system, q is the number of observations in the system S. The fractional part of the parameter is consistent with the fractional part of 
$\delta _{kl}^{GS} $. And code DISB 
$d_{kl}^{GS} $ can be calculated directly.
The ambiguity can be fixed by many methods according to the integer characteristic, such as the LAMBDA (Teunissen, Reference Teunissen1995). The fixed effect of ambiguity can be determined by the ratio (Euler and Schaffrin, Reference Euler and Schaffrin1991), that defined as:
$$\mbox{Ratio}=\frac{\left\| {\hat {N}-{\mathop{N} \limits^{\smile}}^{\prime}} \right\|_{Q_{\hat {N}\hat {N}} }^2 }{\left\| {\hat {N}-{\mathop{N} \limits^{\smile}} } \right\|_{Q_{\hat {N}\hat {N}} }^2 }$$ Where N̂ is the ambiguity floating solution, 
${\mathop{N} \limits^{\smile}}$ indicates the optimal candidate ambiguity group, 
${\mathop{N} \limits^{\smile}}^{\prime}$ is the suboptimal candidate ambiguity group, and 
$Q_{\hat {N}\hat {N}} $ is the variance-covariance (VC) matrix of ambiguity floating solution. 
$\Vert \bullet \Vert _{Q}^{2} $ indicates the Mahalanobis distance that calculated by 
$( \bullet ) ^{\rm T}Q^{-1}( \bullet ) $.
When the ratio value is greater than a certain threshold, the result of the fixed ambiguity is acceptable. If the ambiguity is fixed directly by the calculation result of the Equations (2) and (3), the ratio of the solution result is small due to the influence of the phase DISB. There is a relationship between the phase DISB and the ratio. If a value is subtracted from the inter-system ambiguity result, and then fix the ambiguity. When the value subtracted is the DISB parameter and the ratio with the fixed ambiguity is the largest (Tian et al., Reference Tian, Ge and Neitzel2015, Reference Tian, Ge and Neitzel2017), the principle can be used to calculate the DISB.
Consider the ambiguity calculation result without parameter recombination:
$$\hat {\boldsymbol{N}}=[N_{kl}^{1_g s_1 } ,\ldots ,N_{kl}^{1_g s_{n-1} } ,\bar {\delta }_{kl}^{G{S}_1 } ,\ldots ,\bar {\delta }_{kl}^{G{S}_q } ],\quad \boldsymbol{Q}_{\hat {N}\hat {N}}$$Since the value range of the DISB is −0·5 to 0·5 cycles, the ambiguity floating solution is located on the line segment containing the calculation result of the above formula. Hereinafter, the line, which the line segment locates in, is referred to as a position line, and the endpoints of the line segment are:
$$\left\{ {\begin{array}{@{}l} N_a =[N_{kl}^{1_g s_1 } ,\ldots ,N_{kl}^{1_g s_{n-1} } ,\bar {\delta }_{kl}^{GS_1 } -0.5,\ldots ,\bar {\delta }_{kl}^{GS_q } -0.5] \\ N_b =[N_{kl}^{1_g s_1 } ,\ldots ,N_{kl}^{1_g s_{n-1} } ,\bar {\delta }_{kl}^{GS_1 } +0.5,\ldots ,\bar {\delta }_{kl}^{GS_q } +0.5] \\ \end{array}} \right.$$It can be represented by a parametric equation with the following equation:
$$\boldsymbol{N}(t)=[N_{kl}^{1_g s_1 } ,\ldots ,N_{kl}^{1_g s_{n-1} } ,\bar {\delta }_{kl}^{G{\rm S}_1 } +t,\ldots ,\bar {\delta }_{kl}^{G{\rm S}_q } +t]\quad t\in [-0\cdot 5,0\cdot 5]$$The direction vector of position line is:
$$N_e =[0,\ldots ,0,1,\ldots ,1]$$When the t changes, the point moves on the line segment. If the point moves to a certain position where the ambiguity is fixed to obtain the maximum ratio, and the t can be considered as the phase DISB. However, it is difficult to find the point with the largest ratio value on the line directly.
2.2.1. The point with maximum ratio
The VC matrix of the ambiguity solution can be decomposed by L TDL. The formula of ratio is:
$$\begin{align} \mbox{Ratio}&=\frac{(\hat {N}-{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}}\over{N}} }')^{\rm T}Q_{\hat {N}\hat {N}}^{-1} (\hat{N}-{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}}\over{N}} }')}{(\hat{N}-\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}}\over {N}})^{\rm T}Q_{\hat {N}\hat {N}}^{-1} (\hat{N}-\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}}\over {N}})}=\frac{(\hat{N}-{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}}\over{N}} }')^{\rm T}L^{-1}D^{-1}L^{-\mbox{T}}(\hat{N}-{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}}\over{N}} }')}{(\hat{N}-\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}}\over {N}})^{\rm T}L^{-1}D^{-1}L^{-\mbox{T}}(\hat{N}-\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}}\over {N}})} \notag\\&=\frac{(\hat{N}-{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}}\over{N}} }')^{\rm T}L^{-1}\sqrt {D^{-1}} \sqrt {D^{-1}} L^{-\mbox{T}}(\hat{N}-{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}}\over{N}} }')}{(\hat{N}-\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}}\over {N}})^{\rm T}L^{-1}\sqrt {D^{-1}} \sqrt {D^{-1}} L^{-\mbox{T}}(\hat{N}-\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}}\over {N}})} \notag\\ &=\frac{(\hat{N}-{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}}\over{N}} }')^{\rm T}M^{\rm T}M(\hat{N}-{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}}\over{N}} }')}{(\hat{N}-\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}}\over {N}})^{\rm T}M^{\rm T}M(\hat{N}-\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}}\over {N}})}\end{align}$$
where 
$M=\sqrt {D^{-1}} L^{-{\rm T}}$. The formula shows that the Mahalanobis distance of ambiguity floating solution to the alternative ambiguity group can be transformed into the Euclidean distance in a space, whose base is the column vector of matrix M. Since the matrix M is not an orthogonal matrix necessarily, the space could be tilted. The geometric meaning of the ratio can be expressed as the ratio of the distance between the suboptimal ambiguity candidate group and the optimal ambiguity candidate group in that spatial, and the above conversion is a linear conversion. The direction vector of the position line in this space can be expressed as 
$\textbf{\textit{MN}}_{e} $.
Consider the motion of the point when the t changes. At a certain node position, the point is close to a grid node, and during the subsequent motion, the distance between the motion point and the node gradually decreases. After the distance reaches the minimum value, the motion point will gradually move away from the node until another grid node becomes the nearest one, as shown in Figure 1. It is easy to know that when the point moves to the closest position to the node, the ratio can obtain maximum value, as shown in Figure 2.
Figure 1. Position line in the inclined space and the grid point in the vicinity.
Figure 2. The relationship between the Ratio and the minimum variance and the sub-minimum variance.
According to the geometric principle, when the line connecting the moving point and the grid point is perpendicular to the position line, the grid point is closest to it. At this time, the ratio can be obtained as the optimal candidate ambiguity group. In this case, assume that the optimal ambiguity candidate group in this case is 
${\mathop{N} \limits^{\smile}}$, the ambiguity candidate group needs to be spatially converted first, the candidate group coordinates are 
$M{\mathop{N} \limits^{\smile}}$ in the new space. The formula for calculating the value of t is
$$t=(M{\mathop{N} \limits^{\smile}} -M\hat {N})\cdot MN_e =M({\mathop{N} \limits^{\smile}} -\hat {N})\cdot MN_e $$Where · represents the inner product calculation.
2.2.2. Optimal ambiguity alternative group near the position line
As the position changes during the movement of the motion point on the position line, the optimal ambiguity candidate group will change. All optimal candidate ambiguity groups need to be found and their corresponding t should be calculated separately. Compare the calculation result of the ratio in each case. When the ratio is maximum, the t corresponded is the phase DISB, and the ambiguity fixed result can be treated as the optimal ambiguity result. Since the ambiguity may lose the integer property after converting to another space, it is necessary to find the possible optimal candidate ambiguity group in the original space. In fact, we can use the LAMBDA to find the optimal candidate ambiguity group at any point on the position line. Select some sample points on the position line, and use the LAMBDA to calculate the optimal ambiguity group corresponding to these sample points. When the sample points are dense enough, all possible ambiguity groups can be found. These ambiguity groups correspond to the space grid node. However, too many sample points will result in a huge calculation. In this paper, the following method is used to select the sample point position, to obtain all possible ambiguity groups with as few samples as possible.
In the original space, the range of the position line in each base vector direction is 1. It is difficult to search for the possible ambiguity candidate group due to the influence of correlation. The descending correlation matrix of the ambiguity floating solution variance is Z. After the decorrelation, the direction vector of the position line is:
$$N_{Ze} =Z^{\rm T}N_e$$After the position line transformed into the decorrelation space, the range of variation of the direction in each base vector is N Ze. For the point set on the position line. Obviously, after the optimal ambiguity candidate group changes, the coordinates of the point on a certain base vector will be different unit intervals. The sample point interval is determined by the following formula:
$$\Delta =\frac{1}{\max N_{Ze} }$$Δ is the interval of sample points. After the decorrelation, the formula to calculate the t is as follows:
$$t=M(Z^{\rm T}{\mathop{N} \limits^{\smile}} -Z^{\rm T}\hat {N})\cdot MZ^{\rm T}N_e$$In summary, the method of the phase DISB estimation under the maximum ratio value principle can be divided into the following steps:
1. Calculate the ambiguity floating solution N and its VC-matrix
$Q_{\hat {N}\hat {N}} $, determine the direction vector of the line where the ambiguity floating point solution is located;2. Calculate the decorrelation matrix of the VC-matrix; obtain the interval and the position of the sample points. Get the optimal ambiguity group on each sample point by the LAMBDA;
3. Calculate the t under the optimal ambiguity by the formula;
4. The ambiguity floating solution containing the phase DISB is subtracted from each offset value t and fixed by the LAMBDA to calculate the ratio. The t at the maximum ratio is taken as the result of the final phase DISB.
In addition to the above observation equations and ambiguity fixed strategies, other errors handling strategies and calculation settings in this paper are shown in Table 1. When the ratio is bigger than 2 in the LAMBDA process, the ambiguity is considered to be fixed.
Table 1. Error handling strategy and calculation settings.
3. DATA ANALYSIS
3.1. Data sources and processing methods
Equation 4 represents a tight combination model between the same frequencies and different frequencies. In order to verify the computational performance of the method described in this paper, the data on the GPS L1/Galileo E1, and the data on GPS L1/BDS B1 were analysed respectively. The parameter recombination method and the new method proposed in this paper are used to estimate the DISB. Finally, compared the results of the two methods. The calculated data used short baseline data at Curtin University. The observation time is April 1, 2019, and the relative positions of the stations are shown in Figure 3.
Figure 3. The receivers position in Curtin University.
The different station data were combined to form a baseline combination with the same receivers and different receivers. The baseline information used in the calculation is shown in Table 2.
Table 2. Baselines and the corresponding receiver types in the experiments.
3.2. Analysis of calculation results of DISB between GPS L1 and Galileo E1
There are two types of baselines used to do the data analysis, the first type is the baseline that composed of the same type receivers. And the other is composed of different receivers. In order to test the effect of the method in this paper in calculating the DISB between the systems with the same frequency, we used it to calculate the DISB of the GPS L1 and Galileo E1 data in the above baselines. Then do a comparison with the result of the parameter recombination method. The result series are displayed in the Figure 4. The red point is the result of the parameter combination method and the blue point is the result of the method in the paper in Figure 4. And the (a), (b), (c), (d) are the calculated result of baseline 1, 2, 3, 4.
Figure 4. The DISB results of two methods for 4 baselines. (a) CUT0-CUTB; (b) CUT3-CUBB; (c) CUT0-CUT3; (d) CUTB-CUCC.
Overall, we can get a conclusion that DISB is stable, and if the baseline is composed of the same receivers, it's DISB is close to zero. Additionally, it can be found that code DISB results of the different methods are coincident nearly, which illustrates that the new method cannot influence the calculation of code DISB. Because the method in the paper is based on the regulation in the ambiguity fixed, which just aims to the calculation of phase DISB. However, there is a huge difference in the phase DISB result of two methods in the figure. It can be found that the blue point that representing the new method is much more stable than the red point which representing the parameter recombination method. Indicating that in terms of phase DISB calculation, the results of the method described in this paper can prove the stability of the phase DISB better and get a more reliable result.
There is a convergence process in the blue points in (b) and (c) of Figure 4. It is due to the lack of Galileo data in some epochs. The new method requires several epochs to initialise. In addition, it should be noted that there are some blank parts in (b) and (d) in Figure 4, which is also caused by the lack of Galileo data.
In order to display the result of two methods intuitively, the mean and standard deviation of two methods are calculated and showing in the Table 3.
Table 3. DISB statistical of New Method and Parameter Recombination (PR) Method between GPS L1 and Galileo L1
Obviously, the mean of the DISB between GPS L1 and Galileo E1 is close to 0, which means that DISB can be ignored in the tight combination model between the two frequencies. In the standard deviation part of the code DISB, the new method has no significant improvement effect on the stability. In contrast, the standard deviation of the phase DISB indicates that the method in the paper can improve the stability of phase DISB visibly than the parameter recombination method. The improvement of four baselines are 57·1%, 31·8%, 55·6%, and 43·1%, respectively. Furthermore, as can be seen from the seventh column of Table 3, the method described in this paper can improve the fixed rate of ambiguity to a certain extent. For example, the ambiguity fixed rate increased from 60·36% to 74·27% in the CUTB-CUCC, the promotion is 23·05%. This is because the search domain of the method described in this paper is larger than the parameter recombination method. When phase DISB has a greater influence on the ambiguity fixing, the new method can better overcome the influence of it.
Figure 5 is a histogram of the phase DISB calculation residuals of the two methods. As can be seen from the figure, the residuals of the two methods are within 0·05 cycles. Compared with the calculation residuals of the parameter recombination method, the residuals of the new method are more concentrated, indicating that the calculation results are more stable.
Figure 5. Phase DISB residual of two method for 4 baselines between GPS L1 and Galileo E1. (a) CUT0-CUTB; (b) CUT3-CUBB; (c) CUT0-CUT3; (d) CUTB-CUCC.
3.3. Analysis of calculation results of DISB between GPS L1 and BDS B1
Figure 6 shows the calculation results of the DISB parameters between GPS L1 and BDS B1 frequencies by two methods. It is clear that the phase DISB is stable if regardless of the sudden jump. From the original data analysis, the jumps are due to the data loss of several epochs during observation. It can be inferred from Equation 5 that the result of the phase DISB between different frequencies is related to the ambiguity of the GPS reference satellite in the initial epoch. If there is any loss of observation data, the phase DISB will change as the GPS reference satellite changes after recalculation.
Figure 6. The DISB results of two methods for 4 baselines. (a) CUT0-CUTB; (b) CUT3-CUBB; (c) CUT0-CUT3; (d) CUTB-CUCC.
From the (a) and (b) in Figure 6, code DISB is close to 0 when the baseline is composed of the same receivers even between GPS L1 and BDS B1. But it is far away from 0 if the baseline is composed of different receivers. The value of phase DISB is influenced by the single-differenced ambiguity of the GPS reference satellite and is usually not clear to 0.
In a comparison of the stability part of the two methods' results, the difference of code DISB is small. However, there is a disparity in phase DISB of two methods. The result of the new method is significantly stable than the result of parameter recombination in phase DISB. In addition, as can be seen from (b) and (d) of Figure 6, the new method is more resistant to interference than the parameter recombination method obviously.
Table 4 is a statistical table on DISB between GPS L1 and BDS B1 with two methods. In the CUT0-CUTB and CUTB-CUCC, there are some jumps in the phase DISB series, it does not make sense of calculating the mean value all the time.
Table 4. DISB statistical of New Method and Parameter Recombination (PR) Method between GPS L1 and BDS B1
From the sixth column of the Table 4, the stability improvement of phase DISB is much more obvious than Table 3. The standard deviation of CUT3-CUBB and CUTB-CUCC are reduced from 0·061 cycles and 0·052 cycles to 0·005 and 0·004 cycles, the rate of improvement reach more than 90%. Because in the calculation of the above baselines, there are some gross errors in the result of parameter recombination which make the standard deviation grow. And the new method is more resistant to the errors and the standard deviation of its result is much smaller. For the baseline CUT0-CUT3, the improvement effect of the method described in this paper is not obvious. From the (c) of Figure 6, it can be seen that this is because the calculation results of the two methods are superb.
Figure 7 shows the statistical residuals of the two methods. It can also be found that the histogram of the new method is more concentrated, and it also shows that the stability of the calculation results is stronger. In particular, the concentration of (c) in Figure 7 is not obvious, because the phase DISB time series of CUT0-CUT3 showing a certain time walk characteristics in (c) of Figure 6, so it will cause the concentration of statistical residual influences.
Figure 7. Phase DISB residual of two method for 4 baselines between GPS L1 and BDS B1. (a) CUT0-CUTB; (b) CUT3-CUBB; (c) CUT0-CUT3; (d) CUTB-CUCC.
4. SUMMARY
For the problem that the phase DISB cannot be calculated directly in the tight combination model, a method for the DISB estimation based on the maximum ratio principle was processed in this paper. The method bases that the ambiguity floating solution minus the correct DISB parameter and then calculate by the LAMBDA, the ratio can obtain the maximum value. The ratio calculation formula is used for spatial conversion, and the geometric principle calculation is used to obtain the minimum variance. Finally selected the DISB which makes the ratio maximum as the final result. The observation data in Curtin University were used to verify the method, the following conclusions can be drawn:
From the results of the four sets of baseline calculations, the code DISB and the phase DISB in the tight combination model are stable in one day. The standard deviation of the DISB of the baseline is less than 0·5 m and 0·03 cycles in a day, which composed of the same receivers. About the baselines composed of different receivers, the standard deviation of the code DISB and the phase DISB is smaller than 0·4 m and 0·01 cycles. This result indicates that the DISB can be fixed in advance.
According to the statistics of different baseline calculation results, the ambiguity-fixed rate of the method in the paper is not lower than the parameter recombination method. The promotion can up to 23·05% in some baselines. The reason of the effect is that the ambiguity result of the parameter recombination method is located in the search domain of the method described herein.
Comparing the DISB results between the new method and the parameter recombination. In the code DISB part, the difference between the two methods is small. The improvement of the new method is large only with the poor quality data, the largest promotion in the standard deviation of the code DISB is 0·179 m in the experiment. In the phase DISB part, there is an obvious improvement. All the standard deviation of the phase DISB was decreased with the method in this paper. The maximum rate can up to 92·31% that the standard deviation decreased from 0·052 cycles to 0·004 cycles. This proves the improvement effect in the stability part of phase DISB calculation in this paper.
However, the conclusions above were obtained base on the high-end receivers and short baselines. The application of the method to the low-cost receivers or long baselines should be studied furthermore.
