## 4 GPS/BDS Combined RTK Implementation  
The implementation of the proposed GPS/BDS RTK model is divided into seven modules as Figure 1 shows, and are presented as follows. (1) Single point positioning  
The pseudorange measurements from the rover and the base receivers as well as the broadcast ephemeris decoded by the rover receiver are used as the input to this module. An iterative least square (ILS) method as given in [16] (pp. 122-125) is employed to compute the position of the rover (denoted as \(\mathbf{X}_{\mathrm{r}}\)) and the base (denoted as \(\mathbf{X}_{\mathrm{b}}\)) receivers as well as the elevation angles of all visible satellites. The widely adopted Klobuchar model [17] and Saastamoinen model [18] are used to compensate the ionospheric and tropospheric delays respectively. When processing dual-constellation measurements, we use GPS to generate the single point positions in our implementation, although either constellation or their combination can be selected. It should be pointed out that the BDS time and GPS time have 14 s difference that has to be taken into account when calculating navigation solutions. (2) DD estimation  
This module first generates the non-differenced carrier phase and pseudorange estimation using satellite positions from the ephemeris and the single point rover and base positions from the single point positioning module. The Saastamoinen troposphere model and the mapping function given in [19] are selected to reduce the DD residuals and produce more accurate estimations. Next we single-difference those estimations between the rover and the base, and then double-difference them between the reference satellite(s) and the other satellites. The estimated carrier phases and pseudoranges are denoted as \(\mathbf{\rho}_{\mathrm{DD}}\) and \(\mathbf{\phi}_{\mathrm{DD}}\) respectively. Figure 1: Flowchart of GPS/BDS combined RTK. (3) Measurement processing  
After the measurements have become available as inputs, some processing including detecting cycle slips (e.g., using a loss of lock indicator), handling satellite rising and setting, eliminating satellites below the elevation angle threshold, selecting reference satellite with the largest elevation angle and generating the DD measurements should be done. (4) Pre-filter processing  
This module processes the float SD ambiguity \(\mathbf{N}_{\mathrm{float}}\) (containing both GPS and BDS ambiguities if handling dual-constellation) and the VC matrix \(\mathbf{P}_{k}\) generated from the Kalman filter from the last epoch. If a cycle slip is detected for the \(p\)th satellite, the corresponding column and row of \(\mathbf{P}_{k}\) will be zeroed, \(\mathbf{P}_{k}\)'s corresponding diagonal element will be initialized and the \(p\)th element of \(\mathbf{N}_{\mathrm{float}}\) will be reset to the difference of the SD carrier phase and the SD pseudorange measurements, _i.e._, Equation (7)-Equation (6). The initialization of the filter at the first epoch can be treated as an all-satellite cycle slip case, _i.e._, all \(\mathbf{P}_{k}\)'s diagonal elements and all \(\mathbf{N}_{\mathrm{float}}\)'s elements should be initialized. If one new satellite rises, we should add a row and column to \(\mathbf{P}_{k}\) and initialize the corresponding diagonal element, and \(\mathbf{N}_{\mathrm{float}}\) should be augmented with the corresponding term from Equation (7)-Equation (6) as well. If one satellite that existed at last epoch disappears at the current epoch, then \(\mathbf{P}_{k}\) and \(\mathbf{N}_{\mathrm{float}}\) should be shrunk correspondingly. The measurement matrix \(\mathbf{H}\) is generated in this module based on Equations (20) and (21). At each epoch, we also reset the baseline elements in the state vector to the difference of the single point positions of the rover and the base to avoid dependence on the dynamic features of the base or/and the rover, and reset the corresponding elements in \(\mathbf{P}_{k}\) to a user-defined large number.