## 2 Precise Orbit Positioning Processing Approach  
### Global Navigation Satellite System Observation Model  
The GNSS observation equation of pseudorange \(P\) and carrier phase \(L\) can be expressed as follows:  
\[ \begin{array}{l}L_{r}^{s}=\rho_{r}^{s}+c(dt_{r}-dt^{s})+\lambda(b_{r,s}-b^{s}) +\lambda N_{r}^{s}+T_{r}^{s}-\alpha_{f}I_{r}^{s}+\varepsilon_{r}^{s}\\ P_{r}^{s}=\rho_{r}^{s}+c(dt_{r}-dt^{s})+c(d_{r,s}-d^{s})+T_{r}^{s}+\alpha_{f}I_ {r}^{s}+\varepsilon_{r}^{s}\end{array}  \tag{1}\]  
where \(r\) denotes the receiver, \(s\) the satellite, \(\lambda\) the wavelength of the carrier phase, \(c\) the speed of light in a vacuum, \(\rho_{r}^{s}\) the geometric distance between the satellite position at the signal emission time and receiver position at the signal arrival time, \(dt_{r}\) and \(dt^{s}\) the receiver and satellite clock biases, \(b_{r,s}\) and \(b^{s}\) the receiver- and satellite-dependent UPDs; \(d_{r,s}\) and \(d^{s}\) are the code biases for the receiver and satellite. \(N_{r}^{s}\) is the integer ambiguity, \(T_{r}^{s}\) is the zenith total tropospheric delay, \(I_{r}^{s}\) is the slant ionospheric delay, \(\alpha_{f}\) means the ionospheric delay is a function of the frequency and \(\varepsilon_{r}^{s}\) and \(\varepsilon_{r}^{s}\) are the sum of multipath effect error and measurement noise for the carrier phase and pseudorange observations, respectively. Furthermore, the relativistic delay and phase wind-up errors must be corrected with proper models accordingly [10], although they are not explicitly included in the equations. The ionosphere-free (IF) combination is usually applied to eliminate the first-order ionospheric delays of the code and carrier phase observations. The GNSS code and phase IF combinations are defined as:  
\[ \begin{array}{l}L_{r,IF}^{s}=\frac{f_{1}^{2}}{f_{1}^{2}-f_{2}^{2}}I_{r,i}^{s }-\frac{f_{2}^{2}}{f_{1}^{2}-f_{2}^{2}}I_{r,j}^{s}\\ P_{r,IF}^{s}=\frac{f_{2}^{2}}{f_{1}^{2}-f_{2}^{2}}P_{r,i}^{s}-\frac{f_{2}^{2}}{ f_{2}^{2}-f_{2}^{2}}P_{r,j}^{s}\end{array}  \tag{2}\]  
where \(i\) and \(j\) indicate frequency number and \(f_{i}\) and \(f_{j}\) mean the corresponding frequency. In this study, the L1 and L2 signals are used for GPS and GLONASS, E1 and E5a for Galileo, B1 and B2 for BDS. For simplicity, we will use \(f_{1}\) and \(f_{2}\) to indicate the first and second frequency for all GNSS systems. In POP data processing, satellite and receiver clock offsets are estimated epoch-by-epoch with a network of reference stations with respect to a selected clock as reference. Usually, we can combine the code biases \(d_{r,s}\) and \(d^{s}\) with the clock parameters \(dt^{s}\) and \(dt_{r}\), which turns out to be \(d\overline{t}^{s}\) and \(d\overline{t}_{r}\), respectively. What is more, the code and carrier phases share the same receiver clock, so rewriting Equation (1), we have, for the IF observations,  
\[ \begin{array}{l}L_{r,IF}^{s}=\rho_{r}^{s}+c\left(d\overline{t}_{r}-d \overline{t}^{s}\right)+\lambda_{s,IF}\overline{N}_{r,IF}^{s}+T_{r}^{s}+ \varepsilon_{r,IF}^{s}\\ P_{r,IF}^{s}=\rho_{r}^{s}+c\left(d\overline{t}_{r}-d\overline{t}^{s}\right)+T _{r}^{s}+\varepsilon_{r,IF}^{s}\end{array}  \tag{3}\]  
where through the reformation \where \(r\) and \(r^{\prime}\), \(s\) and \(s^{\prime}\) represent two receivers and two satellites, respectively. The WL ambiguity can be estimated by taking the time average of the Hatch-Melbourne-Wubbena (HMW) [11; 12; 13] combination of the pseudorange and carrier phase observations. Then, the DD WL ambiguity can be fixed according to its probability by rounding to its nearest integer [2; 14]. After WL fixing, the NL ambiguity can be derived with the fixed WL and the float IF ambiguity and fixed in the same way as WL. A DD integer ambiguity can be reconstructed only when both the related WL and NL are fixed, and then as constraints imposed on the normal equations with original UD ambiguities [15].