# Robust Statistics for GNSS Positioning under Harsh Conditions: A Useful Tool?
## 3 Robust Statistics for GNSS Positioning  
The GNSS-based positioning principle consists in solving a geometric problem from the measured ranges to the visible satellites, whose positions are known. Assuming that \(n\geq 4\) satellites are tracked, then the observation model to relate the code pseudoranges to the unknown receiver coordinates is as follows:  
\[R_{i}=\|\mathbf{p}_{i}-\mathbf{p}\|^{2}+\delta t-\delta t_{i}+I_{i}+Tr_{i}+\eta _{i} \tag{24}\]  
where the subscript \(i=\{1,\ldots,n\}\) refers to the \(i\)th satellite, \(R_{i}\) is the observed pseudorange, \(\mathbf{p}_{i}\) and \(\mathbf{p}\) denote the satellite and receiver positions respectively, and \(\delta t\) and \(\delta t_{i}\) are the clock offsets of the receiver and the satellite (in m). In addition, \(I_{i}\) and \(Tr_{i}\) denote the ionospheric and tropospheric corrections and \(\eta_{i}\) gathers the remaining unmodeled errors (e.g., multipath effects, instrumental delays, phase biases, etc.). Solving the system of equations in Equation (24) can be formulated as a regression problem:  
\[\mathbf{y}=h\left(\mathbf{x}\right)+\boldsymbol{\eta} \tag{25}\]  
where \(\mathbf{y}\) is the \(n\)-dimensional observation vector of pseudoranges, \(h\left(\cdot\right)\) is the observation model from Equation (24) and \(\mathbf{x}=\left[\mathbf{p}^{\top},\delta t\right]^{\top}\in\mathbb{R}^{p}\) is the unknown parameter vector. The dimension of the state estimate \(p\) depends on the number of constellations used (three for positioning plus one per each GNSS constellation used). In the context of GNSS, the LS adjustment is the most commonly used method for the estimation of the regression problem of Equation (25). Since GNSS SPP involves a nonlinear observation model, the problem is typically linearized and solved applying an iterative Gauss-Newton method as  
\[\hat{\Delta\mathbf{x}}^{k} =\left(\mathbf{H}^{\top}\mathbf{W}\mathbf{H}\right)^{-1}\mathbf{H}^ {\top}\mathbf{W}\;\left(\mathbf{y}-h\left(\hat{\mathbf{x}}^{k-1}\right)\right) \tag{26}\] \[\hat{\mathbf{x}}^{k} =\hat{\mathbf{x}}^{k-1}+\hat{\Delta\mathbf{x}}^{k} \tag{27}\]  
where \(\mathbf{H}\) is the Jacobian matrix for the observation model, also known as geometry matrix. That linearization is performed around some guess point \(\hat{\mathbf{x}}^{k-1}\) for the \(k\)th iteration of the method, and \(\hat{\Delta\mathbf{x}}\) constitutes the update on the state estimate as in Equation (27). \(\mathbf{W}\) is the weighting matrix for the observations. Classical SPP solutions take \(\mathbf{W}\) as the inverse of the observations covariance matrix \(\mathbf{R}\). Stochastic modelling of pseudorange observations has been a recurrent topic within the GNSS community. A simplification commonly used is to assume that the observations noise is uncorrelated, zero-mean and normally distributed \(\eta_{i}\sim\mathcal{N}\left(0,\sigma_{i}^{2}\right)\)[50]. Thus, the covariance matrix is given by  
\[\mathbf{R}=\mathbf{W}^{-1}=\text{diag}\left(\sigma_{1}^{2},\ldots,\sigma_{n}^ {2}\right) \tag{28}\]  
where \(\sigma_{i}^{2}\) is derived from combining the uncertainty of the different error sources (satellite ephemeris and clock, ionosphere, troposphere, multipath and receiver noise), as in [51; 52] or from error models dependent on the satellite elevation and/or the signal carrier-to-noise density ratio [53; 54; 55]. Algorithm 1 describes the IRLS process for the robust estimation of the GNSS SPP. Notice that WLS (weighted least squares) refers to the iterative Gauss-Newton described in Equations (26) and (27), and MAD is defined in Equation (19). \(N\) and \(\delta\) denote the maximum number of iterations of the iterative Gauss-Newton method and the convergence criteria, respectively. The choice of the influence function and the scale estimate is subject on the robust estimator applied--e.g., for the M-estimator, one might use the Huber function in Equation (10) and the MAD as scale estimate.