## 3 Methodology
### MRER Model  
Based on the previous section, the phase set of multiple satellites selected by different cross-correlation threshold ranges (\(r_{n}\)) can be expressed as:  
\[\begin{cases}x_{a}^{0}=\left[\varphi_{1},\varphi_{2},\ldots,\varphi_{t}\right],(t=1,2,\ldots,m)\\ x=\left[x_{1}^{0},x_{2}^{0},\ldots,x_{a}^{0}\right],(a=1,2,\ldots,32)\end{cases} \tag{8}\]  
where \(x_{a}^{0}\) and \(x\) are the phase set of single satellites and multisatellite combinations, respectively, \(m\) represents the phase length of the satellite, and \(a\) represents the (PRN) number of GNSS satellites.  
The training set samples for building the model were \(\{(x_{i}^{j},y_{i}|i=1,2,\ldots,t_{1};j=1,2,\ldots,\)\(c)\}\)\(\{t_{1}<m,c\leq 32\}\), where \(x_{i}^{j}\) is the multisatellite phase set as the input samples for modeling,is the corresponding SM set as the output samples for modeling, and \(x_{s}^{j}(s=t+1,t+2,\ldots,m)\) is the input sample set for the testing model. Then, a multi-linear regression model was established as follows:  
\[y_{i}=\beta_{0}+\beta_{1}x_{i}^{1}+\beta_{2}x_{i}^{2}+\ldots+\beta_{j}x_{i}^{j}+ v_{i} \tag{9}\]  
where \(i\) represents the sample length of the modeling, \(j\) represents the number of combined satellites for modeling, \(\beta_{0}\), \(\beta_{1}\),..., \(\beta_{j}\) are the regression coefficients of the model, \(v_{i}\) is the residual, including the constant deviation of the regression function and random noise.  
Converting Equation (9) into matrix expression:  
\[\mathbf{Y}=\begin{bmatrix}y_{1}\\ y_{2}\\ \vdots\\ y_{i}\end{bmatrix},\mathbf{X}=\begin{bmatrix}1&\cdots&x_{i}^{1}\\ 1&\cdots&x_{i}^{2}\\ \vdots&\vdots&\vdots\\ 1&\cdots&x_{i}^{j}\end{bmatrix},\mathbf{\beta}=\begin{bmatrix}\beta_{0}\\ \beta_{1}\\ \vdots\\ \beta_{j}\end{bmatrix},\mathbf{V}=\begin{bmatrix}v_{1}\\ v_{2}\\ \vdots\\ v_{i}\end{bmatrix} \tag{10}\]  
where \(Y\) and \(X\) represent the SM matrices and the phase matrices, respectively; \(\mathbf{\beta}\) represents the coefficient matrices that need to be solved; \(\mathbf{V}\) is the residual matrices.  
Therefore, the regression model can be expressed as:  
\[\mathbf{Y}=\mathbf{X}\mathbf{\beta}+\mathbf{V} \tag{11}\]  
A robust estimation can suppress the impact of abnormal values to a certain extent, thus enhancing the number or the quality of the satellites that can be used to monitor SM [63]. In this study, the M estimation was used to replace the traditional least squares estimation. Its basic principle is to use the iterative reweighted least squares (IRLS) method to estimate the regression coefficients [64]. The M estimation criterion can be expressed as:  
\[\sum_{i=1}^{j}\overline{\mathbf{P}}_{ij}\rho(v_{i})=\sum_{i=1}^{j}\overline{\mathbf{P} }_{ij}\rho\big{(}x_{i}\beta_{j}-y_{i}\big{)}=min \tag{12}\]  
where \(\overline{\mathbf{P}}_{ij}\) represents the robust weight function, the selection of which is the key to a robust estimation. In this paper, the IGGIII weight function with \"three-stage\" estimation was adopted, which has strong robustness and can use information more effectively [65; 66].  
Subsequently, according to the IRLS principle, the iterative equation of the robust regression could be obtained as follows:  
\[\mathbf{\beta_{\text{IRLS}}}=\Big{(}\mathbf{X}^{T}\overline{\mathbf{P}}_{ij}\mathbf{X}\Big{)} ^{-1}\mathbf{X}\overline{\mathbf{P}}_{ij}\mathbf{Y} \tag{13}\]  
It can be seen that the main idea of a robust regression model is to assign different weights to different points according to residuals. Then, the weighting algorithm is iterated several times to optimize the weight \(\beta_{\text{IRLS}}\). Therefore, the obtained SM has a certain robustness.