GPS/BDS RTK Positioning Based on Equivalence Principle Using Multiple Reference Stations

Reliable real-time kinematic (RTK) is crucially important for emerging global navigation satellite systems (GNSSs) applications, such as drones and unmanned vehicles. The performance of conventional single baseline RTK (SBRTK) with one reference station degrades greatly in dense, urban environments, due to signal blockage and multipath error. The increasing use of multiple reference stations for kinematic positioning can improve RTK positioning accuracy and availability in urban areas. This paper proposes a new algorithm for multi-baseline RTK (MBRTK) positioning based on the equivalence principle. The advantages of the solution are to keep observation independent and increase the redundancy to estimate the unknown parameters. The equivalent double-differenced (DD) observation equations for multiple reference stations are firstly developed through the equivalent transform. A modified Kalman filter with parameter constraints is proposed, as well as a partial ambiguity resolution (PAR) strategy is developed to determine an ambiguity subset. Finally, the static and kinematic experiments are carried out to validate the proposed algorithm. The results demonstrate that, compared with single global positioning system (GPS) and Beidou navigation system (BDS) RTK positioning, the GPS/BDS positioning for MBRTK can enhance the positioning accuracy with improvement by approximately (45%, 35%, and 27%) and (12%, 6%, and 19%) in the North (N), East (E), and Up (U) components, as well as the availability with improvement by about 33% and 10%, respectively. Moreover, the MBRTK model with two and three reference receivers can significantly increase the redundancy and provide smaller ambiguity dilution of precision (ADOP) values. Compared with the scheme-one and scheme-two for SBRTK, the MBRTK with multiple reference receivers have a positioning accuracy improvement by about (9%, 0%, and 6%) and (9%, 16%, and 16%) in N, E, and U components, as well as the availability improvement by approximately 10%. Therefore, compared with the conventional SBRTK, the MBRTK can enhance the strength of the kinematic positioning model as well as improve the positioning accuracy and availability.


Introduction
Global Navigation Satellite Systems (GNSSs) have been extensively used for scientific and commercial applications in geodesy, geodynamics, transportation, and other industries [1][2][3][4][5]. With the rapid development of GNSSs over the last several decades, the global users will be able to use the multi-constellation and multi-frequency observations to improve the reliability and availability [6][7][8][9][10]. In many emerging GNSS applications, such as precise positioning of drones and unmanned vehicles, the accuracy and availability of real-time kinematic (RTK) are crucially of paramount importance. However, the performance of single-baseline RTK (SBRTK) degrades greatly due to signal blockage and multipath errors.
positioning results of different satellite systems and reference stations are calculated and compared with the SBRTK. This paper is organized as follows. In Section 2, conventional single-differenced (SD) and equivalent double-differenced (DD) observation equations are introduced. Then, a modified Kalman filter with parameter constraints is proposed. The redundancy of MBRTK is given, along with the data processing flow. In Section 3, the static and kinematic experiments are carried out and the impact of satellite systems as well as the reference receivers are discussed. Finally, our conclusions are summarized in Section 4.

SD Observation Equations
As shown in Figure 1, there are four receivers tracking n satellites simultaneously, including three reference receivers r (r = a, b, and c) and one rover receiver u. It is preferable that the precise positions of three reference receivers are obtained in advance. Three reference receivers can be used to enhance the strength of the kinematic positioning model. Without loss of simplicity, the SD observation equations of carrier phase and code observations for single satellite system can be presented as follows: where ( ) ( ) ( ) ru u r       is the SD operator between reference receiver r and rover receiver u. The superscript s (s = j, k, ..., n) denotes the satellite, and satellite j is set as the reference satellite. All satellites can be arranged with the descending order by satellite elevation angle. P and  are the observations of code and carrier phase, respectively, and  denotes the wavelength of the carrier.
 is the geometric distance as a function of the receiver and satellite coordinates. c is the speed of light in vacuum. t represents user receiver clock error. N is the integer ambiguity of carrier phase. T and I are the tropospheric and ionospheric delay, whose difference can be ignored for short baseline.
 and  are the measurement noise of code and carrier phase, respectively. Without loss of simplicity, the SD observation equations of carrier phase and code observations for single satellite system can be presented as follows:  (1) where ∆(·) ru = (·) u − (·) r is the SD operator between reference receiver r and rover receiver u. The superscript s (s = j, k, . . . , n) denotes the satellite, and satellite j is set as the reference satellite. All satellites can be arranged with the descending order by satellite elevation angle. P and ϕ are the observations of code and carrier phase, respectively, and λ denotes the wavelength of the carrier. ρ is the geometric distance as a function of the receiver and satellite coordinates. c is the speed of light in vacuum. t represents user receiver clock error. N is the integer ambiguity of carrier phase. T and I are the tropospheric and ionospheric delay, whose difference can be ignored for short baseline. ς and ε are the measurement noise of code and carrier phase, respectively. By combining SD observation equations of all satellites for single satellite system, the linearized equations in matrix formation can be represented as where x is a column vector of coordinate corrections. N is a column vector of SD ambiguities. t is a column vector of SD receiver clock errors. A and C are the corresponding coefficient matrices, respectively. L is a column vector of the observation minus computation (OMC) terms. v is a vector of SD residual errors. P is a weight matrix of SD observations. The sub-matrices of the above matrices can be calculated by Equation (1) and the specific expressions can be derived as It is assumed that undifferenced observations are independent with elevation-dependent weights and their weight matrix is a diagonal matrix. The different SD observations among multiple receivers tracking the different satellite are not correlative. However, the different SD observations among multiple receivers tracking the identical satellites are correlative, which cannot be ignored. The SD weight matrix can be calculated exactly by the law of variance-covariance propagation. Thus, the weight matrix of Equation (2) is block-diagonal and can be derived as follows where the sub-matrix P s L and P s P are the weight matrix of SD observations for carrier phase and code corresponding to single satellite s.
The specific expression of weight matrix for SD observations of single satellite s is as follows where σ s = a 0 + a 1 e −θ s /10 , σ s is the variance of undifferenced observation, and the symbol θ s is the satellite elevation angle with a 0 = 2 mm, a 1 = 4 mm for carrier phase and a 0 = 20 cm, a 1 = 40 cm for code, the cut-off elevation angle is set to 10 • .

Equivalent DD Observation Equations
In Equation (2), the coordinate corrections and ambiguities are of interest, and the SD receiver clock errors can be eliminated directly from the observation equations based on equivalence principle [31,32]. Thus, the equivalent DD observation equations for single satellite system can be expressed as where E is an identity matrix. R is the transformation matrix of eliminating the clock errors with a rank defect of three. It means that the three SD ambiguities of reference satellite can be set to zero to keep the parameters independent, and the other SD ambiguities can be converted to the equivalent DD ambiguities. Considering the GPS/BDS positioning, the equivalent DD observation equations of MBRTK can be obtained as Remote Sens. 2020, 12, 3178 where the subscripts G and C denote the satellite system.
x is a column vector of coordinate corrections, N G and N C are the column vector of the equivalent DD ambiguities for GPS and BDS systems. As shown in Equation (7), the observations of four receivers can be modelled and processed strictly to obtain the positioning results of the MBRTK, which is not only equivalent to the conventional observation equations of the SBRTK, but also the correlative characteristics of differenced observations are solved clearly. If there are multiple reference and rover receivers tracking multiple satellites, the similar equations can be obtained for RTK positioning. As a special case, when there is only one reference receiver and one rover receiver, the equivalent DD equations can be reduced to the conventional DD equations of the SBRTK.

The Modified Kalman Filter
The Kalman filter is an optimal recursive estimator, which uses a dynamic model and sequential measurements to estimate the unknown state parameters in a minimum variance sense. In this paper, the constant velocity (CV) model [46] are adopted and the position-velocity components of the rover receiver and equivalent DD ambiguities are selected as system states. Let us still take Figure 1 as an example, the system states are expressed as follows: where X au , X bu , and X cu are the position-velocity components of rover receiver, which are estimated corresponding to different reference receivers in Equation (7). In fact, the position-velocity components of rover receiver u for single epoch should be the same values, which can be as a constraint to enhance the strength of RTK positioning. N au , N bu , and N cu are the corresponding equivalent DD ambiguities for three reference receivers, respectively.
Considering the redundant position-velocity components of rover receiver, a constraint equation can be expressed as: Equation (9) can be expanded as follows Thus, based on Equations (7) and (10), the system state and observation equations of Kalman filter with parameter constraints are expressed as follows: where subscript t denotes the epoch, X t is the unknown state vector at epoch t, Φ t,t−1 is the state translation matrix from epoch t − 1 to t, W t is the noise vector of the system state model with zero mean and the covariance matrix Q wk , the setting of covariance matrix Q wk can refer to related literature such as Takasu. [47] and Zhao et al. [17] which is not discussed here. H t is the design matrix and Z t is column vector of the observation minus calculation (OMC) terms at epoch t, V t is a noise vector of observations with zero mean and the covariance matrix Q vk . Moreover, the noise W t and V t are assumed non-correlated.
The Kalman filter solution of the system states with constraint equation [48] is derived as where X t and Q t are the predicted state vector and its corresponding predicted covariance matrix. X t is the corrected value of one-step predicted state vector, which can be improved using constraint equation. J t is the equivalent gain matrix, which can balance the observation and state information,X t and QX t are the state estimation vector and its corresponding covariance matrix at epoch t. With more epoch-data accumulation, the Kalman filter solution will get better and better, and the float ambiguities are more precise in favor to improve the performance of AR.

The Redundancy of MBRTK
The redundancy is an important indicator to illustrate the reliability of the model, which is computed as the number of observations minus the number of unknown parameters for single epoch [49]. In Table 1, we give the number of observations, the number of unknown parameters and the redundancy of MBRTK compared with the SBRTK. Because the Doppler observations are not adopted to estimate the velocity parameters in Equation (7), the velocity components of system states are only predicted in the time update step and are not adjusted in the measurement update step of Kalman filter. Thus, the velocity state parameters are not considered for the calculation of redundancy in Table 1. g and b are the satellite number of single GPS and single BDS. r is the number of reference stations. Because of the equivalence between SBRTK and MBRTK, the condition of MBRTK is the same as the conventional SBRTK, as shown in the fifth column of Table 1. The redundancy is given based on the following assumptions, as (1) the ambiguities are estimated as time-continuous values, (2) the equivalent DD observation equations can be established in intra-system for both single GPS and single BDS, (3) multiple reference receivers can track the same satellites simultaneously. Table 1. The number of observations, unknowns, and redundancy for single epoch where SBRTK = single baseline real-time kinematic, MBRTK = multiple baseline real-time kinematic, GPS = global positioning system, and BDS = Beidou navigation system.

Number of Constraints Condition
Remote Sens. 2020, 12, 3178 7 of 19 As an example, considering the GPS/BDS model for MBRTK, the number of SD observations for both carrier phase and code is 2r(g + b). The number of unknown parameters consists of 3r receiver positions, r(g − 1) ambiguities for GPS, and r(b − 1) ambiguities for BDS. The number of observation redundancy is r(g + b − 1). If there are multiple reference receivers (r ≥ 2) with more data accumulated, the redundancy of MBRTK can increase rapidly, as well as 3(r − 1) constraint equations can be obtained, which can improve the estimable accuracy of unknown state parameters and enhance the strength of MBRTK. Generally, the more the number of reference receivers is, the greater the redundancy improvement is.

Data Processing Strategy
The data processing flow of MBRTK for single epoch is shown in Figure 2, which is divided into two main parts. The first part is data organization and parameter estimation. First, the observations of multiple reference and rover receivers are collected for MBRTK. Second, the SD observations between each reference and rover receivers are used to establish the SD observation equations. Then, the equivalent DD observation equations are obtained through the equivalent transform, and the modified system state equation with parameter constraints are performed to estimate the unknown state parameters. Finally, the float ambiguities and covariance matrices can be derived by the modified Kalman filter.
As an example, considering the GPS/BDS model for MBRTK, the number of SD observations for both carrier phase and code is 2r(g + b). The number of unknown parameters consists of 3r receiver positions, r(g − 1) ambiguities for GPS, and r(b − 1) ambiguities for BDS. The number of observation redundancy is r(g + b − 1). If there are multiple reference receivers (r ≥ 2) with more data accumulated, the redundancy of MBRTK can increase rapidly, as well as 3(r − 1) constraint equations can be obtained, which can improve the estimable accuracy of unknown state parameters and enhance the strength of MBRTK. Generally, the more the number of reference receivers is, the greater the redundancy improvement is.

Data Processing Strategy
The data processing flow of MBRTK for single epoch is shown in Figure 2, which is divided into two main parts. The first part is data organization and parameter estimation. First, the observations of multiple reference and rover receivers are collected for MBRTK. Second, the SD observations between each reference and rover receivers are used to establish the SD observation equations. Then, the equivalent DD observation equations are obtained through the equivalent transform, and the modified system state equation with parameter constraints are performed to estimate the unknown state parameters. Finally, the float ambiguities and covariance matrices can be derived by the modified Kalman filter. The second part of the data processing is AR of MBRTK. A modified PAR method is adopted to determine an ambiguity subset. First, the decorrelated ambiguities rather than raw float ambiguities are chosen, and the first subset is determined by reordering with the ascending estimate precision. Second, the number of valid ambiguities of every single-baseline is calculated. If the number of valid ambiguities is less than 4, only the float solution of MBRTK is available. Third, after the second step, all float ambiguities of MBRTK are used to fix a potential subset by the Agrell, Eriksson, Vardy, Zeger (AEVZ) search algorithm [42]. Fourth, the ratio test and the posterior probability test are calculated. If the subset passes the tests, the current ambiguity subset is accepted with sufficient confidence and the fixed solution of the MBRTK is available. On the contrary, if the last ambiguity with the lowest precision is rejected, the second step is repeated until the ambiguity solution of the MBRTK is determined. The critical criterion of ratio value and posterior probability are set to 2.0 and 0.90, respectively.

Validation and Analysis
In order to test the performance of the proposed MBRTK, we collect the static and kinematic dataset from Curtin University in Australia and Shandong University, Weihai in China, respectively. The static experiment is designed to validate the feasibility and reliability of the proposed MBRTK, and the kinematic experiment is to evaluate its positioning performance for different reference stations. The proposed MBRTK is compared with the conventional SBRTK with PAR method.

Experiment Setup
A one-day dataset on 1 July 2018 of four stations (CUT00, CUTA0, CUTB0, and CUTC0) located at Curtin University in Australia, are collected for the data processing. The positions of four stations are precisely known, and each station is equipped with a TRM59800.00 SCIS antenna, which is connected with Trimble NetR9 receiver (Trimble Navigation Limited, Sunnyvale, CA, USA). Single-frequency code and carrier phase observation of the GPS L1 signal is adopted, and the sampling interval of the observation is 30 s. The satellite cut-off elevation angle is set to 10 • . The distribution of four stations and the sky plot of station CUT00 is depicted in Figure 3. Three schemes are designed to calculate the precise position of CUT00, which are summarized in Table 2. are chosen, and the first subset is determined by reordering with the ascending estimate precision. Second, the number of valid ambiguities of every single-baseline is calculated. If the number of valid ambiguities is less than 4, only the float solution of MBRTK is available. Third, after the second step, all float ambiguities of MBRTK are used to fix a potential subset by the Agrell, Eriksson, Vardy, Zeger (AEVZ) search algorithm [42]. Fourth, the ratio test and the posterior probability test are calculated. If the subset passes the tests, the current ambiguity subset is accepted with sufficient confidence and the fixed solution of the MBRTK is available. On the contrary, if the last ambiguity with the lowest precision is rejected, the second step is repeated until the ambiguity solution of the MBRTK is determined. The critical criterion of ratio value and posterior probability are set to 2.0 and 0.90, respectively.

Validation and Analysis
In order to test the performance of the proposed MBRTK, we collect the static and kinematic dataset from Curtin University in Australia and Shandong University, Weihai in China, respectively. The static experiment is designed to validate the feasibility and reliability of the proposed MBRTK, and the kinematic experiment is to evaluate its positioning performance for different reference stations. The proposed MBRTK is compared with the conventional SBRTK with PAR method.

Experiment Setup
A one-day dataset on 1 July 2018 of four stations (CUT00, CUTA0, CUTB0, and CUTC0) located at Curtin University in Australia, are collected for the data processing. The positions of four stations are precisely known, and each station is equipped with a TRM59800.00 SCIS antenna, which is connected with Trimble NetR9 receiver (Trimble Navigation Limited, Sunnyvale, CA, USA). Singlefrequency code and carrier phase observation of the GPS L1 signal is adopted, and the sampling interval of the observation is 30 s. The satellite cut-off elevation angle is set to 10°. The distribution of four stations and the sky plot of station CUT00 is depicted in Figure 3. Three schemes are designed to calculate the precise position of CUT00, which are summarized in Table 2.    of GPS observed satellites and redundancy is 7.6 and 16.9, respectively. The average RDOP value is 7.3 for the whole day.  Figure 4 depicts the number of observed GPS satellites and redundancy, as well as the relative dilution of precision (RDOP) value of scheme-three during 1 July 2018. As shown, the average number of GPS observed satellites and redundancy is 7.6 and 16.9, respectively. The average RDOP value is 7.3 for the whole day.

Experiment Setup
The kinematic experiment is conducted at Shandong University, Weihai, on November 18, 2018 (from 12:53 to 14:10 UTC). The rover receiver (Trimble NetR9) and antenna (TRM57971.00) are installed on the top of a car, driving on the city roads of Weihai, which is marked as TBR9 as shown in Figure 6. Two reference receivers (Hi-Target VNet) and antenna (HX-AT2300) are located on the rooftop of the Wentian building in Shandong University, while the third reference receiver (South NetS9) and antenna (TRM57971.00) are located on the rooftop of a civil building in Weihai. The three receivers are marked as SDWA, SDWB, and SOTH, respectively. The linear distance of SDWA and SOTH is approximately 6.3 km, and a total number of approximately 4621 epochs of data are collected with 1 s sampling interval and 10 • elevation cut-off angle. The trajectory of the car is shown in Figure 7.   In order to examine the performance of the proposed MBRTK, and compare the impacts of different satellite systems and reference receivers on the positioning results, six different schemes are conducted in this study. The differences between different schemes mainly focus on satellite systems and the number of reference receivers. The description of six positioning schemes is summarized in Table 4. Moreover, the reference coordinates of rover receiver can be calculated by Waypoint Inertial Explorer (IE) software 8.60 (NovAtel, Calgary, Alberta, Canada) for a forward-reverse model with three reference receivers [50], and the positioning errors are defined as the difference between positioning results and the reference coordinates. If the three-dimensional (3D) error of positioning is less than 10 cm, this epoch is defined to be successfully solved for the short baseline, and the availability of kinematic positioning is defined as the percentage of solved epochs of all epochs [38].

Positioning Results of Different Satellite Systems
To evaluate the performance of MBRTK for different satellite systems, the number of observed satellites and redundancy for the scheme-three, scheme-four and scheme-five are shown in Figure 8. In this contribution, G and C are denoted as system identifications for GPS and BDS, and R-is denoted as redundancy. The number of GPS and BDS satellites observed at one epoch over the whole observation segment is more than 6 and 9 mostly, and the total number of GPS/BDS satellites is more than 15. Meanwhile, the average redundancy is 7.6, 13.1, and 26.7 for the three different schemes, respectively. We can see that combined GPS/BDS system can significantly improve the redundancy

Positioning Results of Different Satellite Systems
To evaluate the performance of MBRTK for different satellite systems, the number of observed satellites and redundancy for the scheme-three, scheme-four and scheme-five are shown in Figure 8. In this contribution, G and C are denoted as system identifications for GPS and BDS, and R-is denoted as redundancy. The number of GPS and BDS satellites observed at one epoch over the whole observation segment is more than 6 and 9 mostly, and the total number of GPS/BDS satellites is more than 15. Meanwhile, the average redundancy is 7.6, 13.1, and 26.7 for the three different schemes, respectively. We can see that combined GPS/BDS system can significantly improve the redundancy of kinematic positioning, indicating combined GPS/BDS positioning has a strong ability to test the observations for modeling errors and improve the positioning accuracy. Figure 9 shows the positioning errors of TBR9 after the first convergence for scheme-three, scheme-four, and scheme-five in N, E, and U components, respectively, and the success rate is marked on the graph. The statistics of three schemes are listed in Table 5. The RMS values of three different schemes are (5.5, 2.3, and 6.3) cm, (3.4, 1.6, and 5.7) cm, and (3.0, 1.5, and 4.6) cm in N, E, and U components, while the corresponding standard deviation (STD) values for three schemes are (3.7, 2.3, and 5.5) cm, (2.5, 1.5, and 4.7) cm, and (1.5, 1.4, and 3.9) cm. We can conclude that both the RMS and STD values of combined GPS/BDS system has a better positioning accuracy. Compared with RMS values of scheme-three and scheme-four, the RMS values of scheme-five have a performance improvement of approximately (45%, 35%, and 27%) and (12%, 6%, and 19%) in N, E, and U components, respectively. The STD values of scheme-five also have better consistency than those of scheme-three and scheme-four. Besides, the success rate of scheme-five is 94.4%, which has an availability improvement by approximately 33% and 10% compared with the success rate of scheme-three and scheme-four, respectively.
with RMS values of scheme-three and scheme-four, the RMS values of scheme-five have a performance improvement of approximately (45%, 35%, and 27%) and (12%, 6%, and 19%) in N, E, and U components, respectively. The STD values of scheme-five also have better consistency than those of scheme-three and scheme-four. Besides, the success rate of scheme-five is 94.4%, which has an availability improvement by approximately 33% and 10% compared with the success rate of scheme-three and scheme-four, respectively.   performance improvement of approximately (45%, 35%, and 27%) and (12%, 6%, and 19%) in N, E, and U components, respectively. The STD values of scheme-five also have better consistency than those of scheme-three and scheme-four. Besides, the success rate of scheme-five is 94.4%, which has an availability improvement by approximately 33% and 10% compared with the success rate of scheme-three and scheme-four, respectively.

Positioning Results of Different Reference Receivers
To evaluate the performance of MBRTK for different number of reference receivers, the redundancy for the scheme-one, scheme-two, scheme-five and scheme-six are shown in Figure 10. The average redundancy for scheme-one and scheme-two is 9.4 and 9.3, which are almost the same because of the same processing strategy with the traditional SBRTK model, while the average redundancy for scheme-five and scheme-six with MBRTK is 26.7 and 40.0. The more the reference receivers are, the stronger the reliability of MBRTK model is. Meanwhile, the ADOP values [51,52] quantify the priori precision and geometry of ambiguities, which can be computed and illustrated in Figure 11 for the four schemes. We can see that, the ADOP values of all the four schemes can remain below 0.30 and 0.15 cycles for 99% and 90% of the observation segment, respectively. The ADOP values of scheme-five and scheme-six is the smaller than those of scheme-one and scheme-two, and scheme-six with three reference receivers has the optimal ADOP values. This is, thus, a promising indication that faster successful AR is possible for MBRTK model using multiple reference receivers. below 0.30 and 0.15 cycles for 99% and 90% of the observation segment, respectively. The ADOP values of scheme-five and scheme-six is the smaller than those of scheme-one and scheme-two, and scheme-six with three reference receivers has the optimal ADOP values. This is, thus, a promising indication that faster successful AR is possible for MBRTK model using multiple reference receivers. Figure 12 shows the positioning errors of TBR9 after the first convergence for scheme-one, scheme-two, and scheme-six in N, E, and U components, respectively. Considering the positioning errors of scheme-five discussed in the previous section, the corresponding statistics of four schemes with different reference receivers are listed in Table 6 Figure 12 shows the positioning errors of TBR9 after the first convergence for scheme-one, scheme-two, and scheme-six in N, E, and U components, respectively. Considering the positioning errors of scheme-five discussed in the previous section, the corresponding statistics of four schemes with different reference receivers are listed in Table 6. The RMS values of four different schemes are     From the analysis of the RMS and STD values, we can conclude that it is no statistically significant difference in the positioning errors for scheme-five and scheme-six, which can achieve a similar accuracy of 3.0 cm, 1.5 cm, and 4.5 cm level in N, E, and U components. Compared with the RMS values of scheme-one and scheme-two with SBRTK model, the RMS values of scheme-five and scheme-six with MBRTK model have a performance improvement by approximately (9%, 0%, and 6%) and (9%, 16%, and 16%) in N, E, and U components, respectively. Meanwhile, the STD values of scheme-five and scheme-six illustrate a better consistency than those of scheme-one and scheme-two. Besides, the success rate of four different schemes is 84.5%, 86.0%, 94.4% and 96.8%, respectively. Compared with scheme-one and scheme-two, the scheme-five and scheme-six have a similar From the analysis of the RMS and STD values, we can conclude that it is no statistically significant difference in the positioning errors for scheme-five and scheme-six, which can achieve a similar accuracy of 3.0 cm, 1.5 cm, and 4.5 cm level in N, E, and U components. Compared with the RMS values of scheme-one and scheme-two with SBRTK model, the RMS values of scheme-five and scheme-six with MBRTK model have a performance improvement by approximately (9%, 0%, and 6%) and (9%, 16%, and 16%) in N, E, and U components, respectively. Meanwhile, the STD values of scheme-five and scheme-six illustrate a better consistency than those of scheme-one and scheme-two. Besides, the success rate of four different schemes is 84.5%, 86.0%, 94.4% and 96.8%, respectively. Compared with scheme-one and scheme-two, the scheme-five and scheme-six have a similar availability improvement by approximately 10% in the success rate. However, it should be noted that scheme-six with three reference receivers does not present better performance than scheme-five with two reference receivers in both positioning accuracy and availability. The possible reason is that the reference receiver SDWB is close to the SDWA, which cannot further enhance the strength of the control network in the kinematic positioning. Figure 13 presents L1/B1 carrier residuals for four schemes with different reference receivers. The DD residuals can be derived with the fixed coordinate parameters for MBRTK to compare with the residuals of SBRTK. It can be seen that the residuals of four schemes are randomly distributed along the zero Y-axis, and the mean values of the residuals are approximately 0 cm. The RMS and STD values of the residuals for scheme-one are worse than that of scheme-two. The possible reason is that there is a mountain behind the SDWA station, which leads to a worse satellite geometry than that of SOTH station. Although the residuals of scheme-six are slightly worse than that of scheme-five, the RMS and STD values of both schemes with multiple reference receivers have a good consistency, which are obviously smaller than that derived from scheme-one and scheme-two.
along the zero Y-axis, and the mean values of the residuals are approximately 0 cm. The RMS and STD values of the residuals for scheme-one are worse than that of scheme-two. The possible reason is that there is a mountain behind the SDWA station, which leads to a worse satellite geometry than that of SOTH station. Although the residuals of scheme-six are slightly worse than that of schemefive, the RMS and STD values of both schemes with multiple reference receivers have a good consistency, which are obviously smaller than that derived from scheme-one and scheme-two.

Discussion
The results of Section 3 indicate that the reference coordinates of rover receiver (TBR9) are calculated by Waypoint IE software 8.60 for a forward-reverse model with three reference receivers. However, the forward Kalman filter is applied to obtain the positioning results in MBRTK algorithm. The differences of model in the data processing would lead to a systematic positioning error, as shown in Table 6.
In addition, the positioning results of the forward-reverse model of IE software 8.60 still have the positioning error compared to the true position of rover receiver. Thus, to study the above

Discussion
The results of Section 3 indicate that the reference coordinates of rover receiver (TBR9) are calculated by Waypoint IE software 8.60 for a forward-reverse model with three reference receivers. However, the forward Kalman filter is applied to obtain the positioning results in MBRTK algorithm. The differences of model in the data processing would lead to a systematic positioning error, as shown in Table 6.
In addition, the positioning results of the forward-reverse model of IE software 8.60 still have the positioning error compared to the true position of rover receiver. Thus, to study the above differences, three additional schemes with different single reference receiver are conducted using the forward-reverse model by IE software 8.60. The statistics are listed in Table 7. As can be seen, both the RMS and STD values of SDWA-TBR9 for IE software 8.60 can achieve the mm level and the positioning results are approximately equal to the reference coordinates of rover Remote Sens. 2020, 12, 3178 16 of 19 receiver by IE software 8.60. However, for the RMS and STD values of SDWB-TBR9 and SOTH-TBR9 using IE software 8.60, there is an obvious systematic error with cm level. That is to say, the reference coordinates by IE software 8.60 with three reference receivers are probably weighted average values based on the positioning results of SDWA-TBR9, SDWB-TBR9, and SOTH-TBR9, and the weight of SDWA-TBR9 is largest among the three single baselines. In other words, the reference coordinates of TBR9 are only relatively precise, but not absolutely accurate. It is worth noting that although the scheme-five and scheme-six still have a systematic error with cm level, the model of MBRTK is rigorous and reliable, which also have a higher success rate than those of three additional schemes by IE software 8.60.

Conclusions
In this paper, we present a new algorithm of MBRTK based on the equivalence principle for precise kinematic positioning. First, equivalent DD observation equations using SD observations are obtained through the equivalent transform to eliminate the receiver clock errors. Second, considering the redundancy of state information, a modified Kalman filter with parameter constraints is proposed to estimate the position-velocity states and equivalent DD ambiguities. The most outstanding feature of the MBRTK is that it has a rigorous model to solve the correlative characteristics of differenced observations for multiple reference receivers, as well as high positioning accuracy and availability could be achieved.
From the experiment and analysis, the following conclusions can be obtained: (1) For the static experiment, the MBRTK has better positioning accuracy than SBRTK. The MBRTK results with two and three reference receivers have the improvement of positioning accuracy by approximately (14%, 13%, and 14%), and (21%, 17%, and 19%) in N, E, and U components, respectively, validating the feasibility and reliability of the proposed MBRTK. (2) For the kinematic experiment, combined GPS/BDS positioning for MBRTK increases the number of available satellites and redundancy of the positioning model. Compared with single GPS and single BDS positioning, the combined GPS/BDS positioning has the accuracy improvement by approximately (45%, 35%, and 27%) and (12%, 6%, and 19%) in N, E, and U components and the availability improvement by approximately 33% and 10%, respectively. (3) For the analysis of reference receivers, the MBRTK model with multiple reference receivers can significantly increase the redundancy and provide the smaller ADOP values in favor to improve the performance of AR in comparison to the SBRTK model. The positioning results with two and three reference receivers can provide the similar accuracy of 3.0 cm, 1.5 cm, and 4.5 cm level in N, E, and U components, which have the accuracy improvement by approximately (9%, 0%, and 6%) and (9%, 16%, and 16%) compared with two SBRTK results, respectively. Meanwhile, the MBRTK results with two and three reference receivers have the availability improvement by approximately 10% in the success rate. (4) Due to the different models adopted by IE software 8.60 and MBRTK, there exists systematic errors between them. Compared with the positioning results of IE software 8.60 with one reference receiver, the MBRTK can achieve the stable accuracy with the cm-level, and have a rigorous and reliable model to keep a high availability of kinematic positioning.
Our future work will extend this algorithm to multi-frequency and multi-GNSS positioning considering the computational efficiency, and focus the distribution of multiple reference receivers to improve the precision, availability, and reliability of positioning results.