## 2 Methods  
In this section, the multisystem integrated SPP model taking ISB as a priori constraints is described. Only the code observations on E1 frequency for Galileo, on B1 frequency for BDS, on G1 frequency for GLONASS, and on L1 frequency for GPS are adopted here, since most of the SPP users are still using the low-cost single-frequency receivers. The single-frequency code observations from a GNSS satellite can be modeled as:  
\[P=\rho+c\mathrm{d}t_{r}-c\mathrm{d}t+I+T+d_{r}+d \tag{1}\]  
where \(P\) is the measured pseudorange, \(\rho\) is the geometric range between the phase centers of the satellite and receiver antennas, \(c\mathrm{d}t_{r}\) and \(c\mathrm{d}t\) are the physical clock errors of the receiver and the satellite, respectively, \(I\) is the slant ionospheric delay, \(T\) is the slant tropospheric delay, and \(d_{r}\) and \(d\) are the stable code hardware delays at the receiver and the satellite, respectively. The alignment of time reference and coordinate reference is one of key issues when performing multi-GNSS integrated data processing. The coordinate transformations can be omitted for the SPP processing with a meter-level accuracy, as the difference of coordinate references among different satellite systems (i.e., GTRF for Galileo, CGCS2000 for BDS, PZ90.11 for GLONASS, and WGS-84 for GPS) is confined to several centimeters [16]. By contrast, the difference of the time scales employed by the four GNSSs (i.e., Galileo System Time for Galileo, BDS Time for BDS, GLONASS System Time for GLONASS, and GPS Time for GPS) must be properly considered in the multi-GNSS combined SPP, in view that the effect of the time reference is significant. To mitigate the effect of distinct time scales in the multi-GNSS combined SPP, one receiver clock offset parameter can be estimated for each satellite system. Alternatively, by taking the receiver clock offset parameter of a selected satellite system as the reference, for example, GPS, the additional ISB parameters can be introduced into the code observation equations for other satellite systems. In addition to the distinct time scales, the ISB parameters can also absorb the effect of different receiver-specific code hardware delays among different GNSSs. Due to the inconsistent code hardware delays between BDS-2 and BDS-3 satellites [10; 11], an ISB parameter for BDS-2 and an ISB parameter for BDS-3 are estimated. After correcting the ionospheric delay and tropospheric delay with a priori model, fixing the satellite position with broadcast satellite orbits, and correcting the satellite clock offset with broadcast satellite clocks and time group delays (TGDs), the linearized observation model for four-system integrated SPP can be expressed as:  
\[\left\{ \begin{array}{l}p^{G}=\mu^{G}\cdot X+c\mathrm{d}t_{r,est}\\ p^{R}=\mu^{R}\cdot X+c\mathrm{d}t_{r,est}+\mathrm{ ISB}^{R,G}\\ p^{E}=\mu^{E}\cdot X+c\mathrm{d}t_{r,est}+\mathrm{ ISB}^{E,G}\\ p^{C2}=\mu^{C2}\cdot X+c\mathrm{d}t_{r,est}+\mathrm{ ISB}^{C2,G}\\ p^{C3}=\mu^{C3}\cdot X+c\mathrm{d}t_{r,est}+\mathrm{ ISB}^{C3,G}\\ \end{array} \right. \tag{2}\]  
where the superscripts \(G\), \(R\), \(E\), C2, and C3 refer to GPS, GLONASS, Galileo, BDS-2, and BDS-3 satellites, respectively, \(p\) is the observed-minus-computed (OMC) code observable, \(\mu\) is the unit vector in line-of-sight direction, \(X\) is the receiver coordinates in three dimensions (corrections referring to the approximate values), \(c\mathrm{d}t_{r,est}\) is the receiver clock estimate grouped with the receiver-specific code hardware delays of GPS, and \(\mathrm{ ISB}^{R,G}\)ISB\({}^{E,G}\), ISB\({}^{\text{C2},G}\), and ISB\({}^{\text{C3},G}\) denote the ISB parameter between GLONASS and GPS, between Galileo and GPS, between BDS-2 and GPS, and between BDS-3 and GPS, respectively. Actually, a priority list should be made for the selection of reference satellite system for the ISB estimation.