## 3 Robust Extended Kalman Filter
### RAIM Algorithm  
In this paper, RAIM algorithm is used to detect positioning errors exceeding the alert limit. One of RAIM's outputs is the HPL, that is defined as the radius of a circle in the horizontal plane, centered at the true position. The Horizontal Alert Limit (HAL) is the maximum allowed HPL, which means that the estimated positions are accepted when \(HPL\leq HAL\). Given the linear Equation (27), the estimated state vector \(\hat{X}\) is determined by IRWLS as explained in Section 3.3:  
\[\hat{X}=\left(G^{T}WG\right)^{-1}G^{T}Wy \tag{41}\]  
where \(W\) is the \(W\) matrix last obtained matrix in Equation (40). The state error is calculated as the difference between the estimated state vector \(\hat{X}\) and the true state vector \(X\) :  
\[ \begin{split}\delta X=\hat{X}-X&=\left(G^{T}WG \right)^{-1}G^{T}Wy-X=\left(G^{T}WG\right)^{-1}G^{T}W(GX+e)-X\\ &=\left(G^{T}WG\right)^{-1}G^{T}We\end{split}  \tag{42}\]  
The estimated residual \(\hat{e}\) is determined by computing the difference between the measured observation \(y\) and the estimated observation (\(G\hat{X}\)):  
\[ \begin{split}\hat{e}=y-G\hat{X}=y&-G\left(G^{T}WG \right)^{-1}G^{T}Wy=\left(I-G\left(G^{T}WG\right)^{-1}G^{T}W\right)y=Sy\\ &=\left(I-G\left(G^{T}WG\right)^{-1}G^{T}W\right)(GX+e)\\ &=GX+e-G\left(G^{T}WG\right)^{-1}G^{T}WGX-G\left(G^{T}WG\right)^{ -1}G^{T}We\\ &=GX+e-GX-G\left(G^{T}WG\right)^{-1}G^{T}We=\left(I-G\left(G^{T}WG \right)^{-1}G^{T}W\right)e=Se\end{split}  \tag{43}\]  
where:  
\[S=I-G\left(G^{T}WG\right)^{-1}G^{T}W \tag{44}\]Note that \(S\) is a projector matrix and thus idempotent: \(S^{2}=S\). The a priori standard deviation of the pseudo-range is calculated as follows:  
\[\sigma_{0}=\sqrt{\hat{e}^{T}W\hat{e}}/\sqrt{n-p} \tag{45}\]  
where \(n\) is the number of observations and \(p\) is the number of parameters in the unknown state \(X\). The position error in the state error \((\delta X)\) of Equation (42) is defined in the (x, y, z) coordinate system. In practice, a local (e, n, u: east, north, up) coordinate system is more suitable. To transform from (x, y, z) coordinates to (e, n, u) coordinates, the \(F^{T}\) orthogonal transformation matrix is required [23]:  
\[F^{T}=\left[ \begin{array}{ccc}-\sin\lambda&\cos\lambda&0\\ -\sin\varphi\cos\lambda&-\sin\varphi\sin\lambda&\cos\varphi\\ \cos\varphi\cos\lambda&\cos\varphi\sin\lambda&\sin\varphi\end{array} \right] \tag{46}\]  
where \(\varphi\) is the longitude and \(\lambda\) the latitude. Here, only the parameters of coordinates are considered. For example, the state vector \(X\) is defined following Equation (10) with the first, third and fifth positions in relation to the three (x, y, z) coordinates. Therefore, the new \(\tilde{\mathcal{G}}_{0}\) matrix is obtained as \(\tilde{\mathcal{G}}_{0}=G(;1:2:5)\) and the position errors as \(\delta\tilde{X}_{0}=\delta X(1:2:5)\). The \(M\) matrix related to the errors in easting (\(dE\)), northing (\(dN\)) and upping (\(dU\)) is calculated as follows:  
\[\delta X_{ENU}=\left[ \begin{array}{c}\mathrm{d}E\\ \mathrm{d}N\\ \mathrm{d}U\end{array} \right]=F^{T}\delta\tilde{X}_{0}=F^{T}\left(\mathcal{G }_{0}^{T}\mathcal{G}_{0}\right)^{-1}\tilde{G}_{0}^{T}e=Me \tag{47}\]  
Since the horizontal error related to easting and northing is considered, a matrix \(\tilde{M}_{0}\) can be defined as \(\tilde{M}_{0}=M(1:2,:)\). In reality, none of the elements in \(e\) is equal to zero. However, when serious observation measure error occurs, some elements in \(e\) may be extremely larger than the others. Therefore, we only consider the \(ith\) observation that has a failure of magnitude \(\beta\), while we suppose that the other elements are equal to zero. The error vector can thus be re-written as:  
\[e=\left[ \begin{array}{cccc}0&\cdots&\beta&\ldots&0\end{array}