## Methodology  
### PPP Model and Data Processing Strategy  
After correcting the satellite orbits and clocks using the precise products, the PPP model for GPS, GLONASS, and BDS can be written as follows:  
\[\left\{ \begin{array}{l}p_{r,j}^{G}=-\mathbf{e}_{r}^{G}\cdot\mathbf{r}_{r}+ \bar{\mathbf{t}}_{r}+\gamma_{j,G}\overline{I}_{r,1}^{G}+T_{r}^{G}+M_{r,j}^{G}+ \xi_{r,j}^{G}\\ p_{r,j}^{Rk}=-\mathbf{e}_{r}^{Rk}\cdot\mathbf{r}_{r}+\bar{\mathbf{t}}_{r}+ \delta_{r,Rk}+\gamma_{j,Rk}^{R}\bar{I}_{r,1}^{Rk}+T_{r}^{R}+M_{r,j}^{Rk}+\xi_ {r,j}^{Rk}\\ p_{r,j}^{C}=-\mathbf{e}_{r}^{C}\cdot\mathbf{r}_{r}+\bar{\mathbf{t}}_{r}+\delta _{r,C}+\gamma_{j,C}\overline{I}_{r,1}^{C}+T_{r}^{C}+M_{r,j}^{C}+\xi_{r,j}^{C} \end{array} \right., \tag{1}\]  
\[\left\{ \begin{array}{l}l_{r,j}^{G}=-\mathbf{e}_{r}^{G}\cdot\mathbf{r}_{r}+ \bar{\mathbf{t}}_{r}-\gamma_{j,G}\overline{I}_{r,1}^{G}+T_{r}^{G}+\lambda_{j}^ {C}\overline{N}_{r,j}^{G}+m_{r,j}^{G}+\varepsilon_{r,j}^{G}\\ l_{r,j}^{Rk}=-\mathbf{e}_{r}^{Rk}\cdot\mathbf{r}_{r}+\bar{\mathbf{t}}_{r}+\delta _{r,Rk}-\gamma_{j,Rk}\overline{I}_{r,1}^{R}+T_{r}^{R}+\lambda_{j,Rk}^{R} \overline{N}_{r,j}^{R}+m_{r,j}^{Rk}+\varepsilon_{r,j}^{R}\\ l_{r,j}^{C}=-\mathbf{e}_{r}^{C}\cdot\mathbf{r}_{r}+\bar{\mathbf{t}}_{r}+\delta _{r,C}-\gamma_{j,C}\overline{I}_{r,1}^{C}+T_{r}^{C}+\lambda_{j}^{C}\overline{ N}_{r,j}^{C}+m_{r,j}^{C}+\varepsilon_{r,j}^{C}\end{array} \right., \tag{2}\]  
where \(p_{r,j}^{s}\) and \(l_{r,j}^{s}\) are the \"observed minus computed\" pseudo-range and phase observations from receiver \(r\) to satellite \(s\) on frequency \(j\) (\(j=1,2\)); the superscripts G, R, and C represent the GPS, GLONASS, and BDS systems, respectively; \(R_{k}\) denotes the frequency factor of the GLONASS satellite; \(\mathbf{e}_{r}^{s}\) is the line-of-sight unit vector from receiver to satellite; \(\mathbf{r}_{r}\) denotes the vector of the receiver position increments with respect to the a priori position used for linearization; \(\bar{\mathbf{t}}_{r}\) is the receiver clock offset; \(\delta_{r,C}\) and \(\delta_{r,Rk}\) denote the inter-system bias (ISB) with respect to GPS and the inter-frequency bias (IFB) for GLONASS, respectively; and \(\overline{N}_{r,j}^{s}\) refers to the float ambiguity containing code hardware and phase delay, while \(\lambda_{j}^{s}\) is its corresponding wavelength. The ionospheric delay at different frequencies can be expressed as \(\overline{I}_{r,j}=\gamma_{j,s}\overline{I}_{r,1},\gamma_{j,s}=f_{1}^{2}/f_ {j}^{2}\); \(T_{r}^{s}\) denotes the slant tropospheric delay; \(M_{r,j}^{s}\) and \(m_{r,j}^{s}\) represent the multipaths in pseudo-range and carrier-phase observations; and \(\xi_{r,j}^{s}\) and \(\varepsilon_{r,j}^{s}\) are the measurement noises. Other error items, such as the phase wind-up, relativity effects, Earth rotation, and tidal loading, are corrected by applying models described in Kouba [27]. It should be noted that, for GPS and GLONASS, the phase center offsets and variations (PCO and PCV) at both satellite and receiver are obtained from the International GNSS Service (IGS) antenna file. On the other hand, for BDS, the PCO and PCV at satellite are available from the IGS antenna file, but replaced by GPS at the receiver [28]. The data processing information is listed in Table 1 in detail. The receiver clock offset is treated as white noise, and estimated epoch by epoch. Together with other parameters, the receiver positions are estimated as daily solution using the least-squares estimator. For the SF method, the satellite orbital repeat times for GPS and BDS satellites are computed individually using Keplerian orbital elements from broadcast ephemerides (19; 31). On the other hand, for GLONASS, its broadcast ephemeris is presented by positions and velocities, and thus we used the aspect repeat time instead (32; 33). The phase residuals over n days before the day of interest are shifted by n times the orbital repeat time and then stacked for each station-satellite pair.