## 2 Methodology
### Observation Model of Multi-GNSS PPP  
For a satellite s observed by a receiver r, the basic equations for raw GNSS pseudorange and carrier phase observations are expressed as follows:  
\[ \begin{split}& p^{s}_{r,j}=\mu^{s}_{r}\Delta x_{r}+t_{r}-t^{s}+d^{ s}_{r,j}-d^{s}_{j}+\beta^{s}_{r,j}I_{z,r}+m^{s}_{r}Z_{r}+\varepsilon^{s}_{r,j}\\ & l^{s}_{r,j}=\mu^{s}_{r}\Delta x_{r}+t_{r}-t^{s}-\beta^{s}_{r,j} I_{z,r}+m^{s}_{r}Z_{r}+\lambda^{s}_{j}\Big{(}B^{s}_{r,j}+b^{s}_{r,j}-b^{s}_{j} \Big{)}+\xi^{s}_{r,j}\end{split}  \tag{1}\]  
where \(p^{s}_{r,j}\) and \(l^{s}_{r,j}\) are the observed minus computed (OMC) pseudorange and phase observations on frequency band j in meters; \(\mu^{s}_{r}\) is the unit row vector from the receiver to the satellite after linearization; \(\Delta x_{r}\) is the column vector of the receiver's position increments relative to the a priori position; \(t_{r}\) and \(t^{s}\) are the clock offsets for the receiver and satellite, respectively; \(d^{s}_{r,j}\) and \(d^{s}_{j}\) are the frequency-dependent uncalibrated code delays (UCDs) for the receiver and satellite; \(I_{z,r}\) is the zenith total ionospheric content and \(\beta^{s}_{r,j}=\gamma^{s}_{T}\cdot 40.3/f^{2}_{j}\) is a frequency-dependent scalar factor, in which \(\gamma^{s}_{r}\) is the ionospheric mapping function and \(f_{j}\) is the frequency [30]; \(Z_{r}\) is the zenith tropospheric delay and \(m^{s}_{r}\) is the corresponding mapping function; \(B^{s}_{r,j}\) is the integer phase ambiguity in unit of cycles and \(\lambda^{s}_{j}\) is the corresponding wavelength; \(b^{s}_{r,j}\) and \(b^{s}_{j}\) are the frequency-dependent UPDs for the receiver and satellite; and \(\varepsilon^{s}_{r,j}\) and \(\xi^{s}_{r,j}\) are the measurement noise of pseudorange and carrier phase observations. The ionospheric-free pseudorange and carrier phase observables are usually utilized to remove the first-order ionospheric delay errors [31]. The satellite-dependent UCDs can be also eliminated from the ionospheric-free observation model as it is involved in the precise satellite clock products, which are generated with the ionospheric-free undifferenced observations [32]. Because of different signal frequencies and signal structures, the intersystem biases (ISBs) and interfrequency biases (IFBs) should be considered for multi-GNSS combined PPP [33]. However, the introduction of the GLONASS IFB parameter will increase the complexity of the solution. Considering that it does not influence the validation of this proposed robust approach, the GLONASS IFB parameter is approximately treated as one ISB parameter. If GPS system is chosen as the reference system, the ionospheric-free model of the quad-system combination may be expressed as  
\[ \begin{split}& p^{s,G}_{r,IF}=\mu^{s,G}_{r}\Delta x_{r}+\bar{t}_{r}+m^{s,G}_{r}Z _{r}+\varepsilon^{s,G}_{r,IF}\\ & l^{s,G}_{r,IF}=\mu^{s,G}_{r}\Delta x_{r}+\bar{t}_{r}+\bar{B}^{ s,G}_{r,IF}+m^{s,G}_{r}Z_{r}+\varepsilon^{s,G}_{r,IF}\end{split}  \tag{2}\]  
\[ \begin{split}& p^{s,SYS}_{r,IF}=\mu^{s,SYS}_{r}\Delta x_{r}+\bar{t}_{r}+ ISB^{SYS-G}_{r}+m^{s,SYS}_{r,IF}Z_{r}+\varepsilon^{s,SYS}_{r,IF}\\ & l^{s,SYS}_{r,IF}=\mu^{s,SYS}_{r}\Delta x_{r}+\bar{t}_{r}+ ISB^{SYS-G}_{r}+\overline{B}^{s,SYS}_{r,IF}+m^{s,SYS}_{r}Z_{r}+\varepsilon^{s,SYS}_{r,IF}\end{split}  \tag{3}\]  
where \"\(SYS\)\" represents the GLONASS (R), BDS (C), and Galileo (C) systems, respectively, and the symbol \"\(G\)\" refers to the GPS system; the receiver clock can be expressed as  
\[\bar{t}_{r}=t^{G}_{r}+d^{G}_{r,IF} \tag{4}\]  
and the biases are defined as  
\[ISB^{SYS-G}_{r}=d^{SYS}_{r,IF}-d^{G}_{r,IF} \tag{5}\]where \(d^{G}_{r,IF}\) and \(d^{SYS}_{r,IF}\) denote the ionospheric-free receiver UCDs of the GPS reference system and other systems. Considering the satellite clock products contain the satellite-dependent UCDs, the ambiguity can be expressed as  
\[