## 6 Proposed Multiconstellation Positioning Algorithm
### \(X_{0}\) Estimation (Stage 1)  
1. Select \(P\) subsamples of size \(n_{sub}\) from all \(n+m\) observations. The number of the subsamples is, \[\overline{N}=\left(\begin{array}{c}n+m\\ n_{sub}\end{array}\right),\hskip 28.452756 pt\#(P)=n_{sub}\] (23)  
\(n_{sub}\) represents the number of observations in one subsample. The value \(n_{sub}\) should always meet the criteria of \(g_{n}<n_{sub}<n+m\). \(g_{n}\) is the minimal number of observations required for positioning. \(n_{sub}\) is selected so that at least one combination contains all observations with high \(CNo\). However, Equation (23) is not directly applicable to a multiconstellation positioning system. For a multiconstellation positioning system, the number of unknowns are \(3+\Re\), here \(\Re\) is the number of constellations used in positioning (clock correction parameters). In order to attain the benefits of the multiconstellation positioning, subsamples must contain at least one observations from each constellation, for example a two constellations based positioning system with five observations independently, maximum four observations can be selected from either constellation to fulfill the minimum requirement as minimum five observations are required for positioning.
2. Let \(p_{k}\) is a subsample belonging to \(P\) where \(P=\big{\{}p_{1},p_{2},\ldots,p_{\overline{N}}\big{\}}\). We compute the initial estimate \(X_{p_{k}}^{0}\) with the WLSE method, \[X_{p_{k}}^{0}=\Big{(}\mathbf{H}_{p_{k}}^{T}\Sigma_{p_{k}}^{0}\mathbf{H}_{p_{k }}\Big{)}^{-1}\mathbf{H}_{p_{k}}^{T}\Sigma_{p_{k}}^{0}\rho_{p_{k}}\] (24) Here \(\mathbf{H}_{p_{k}}\), \(\rho_{p_{k}}\) and \(\Sigma_{p_{k}}^{0}\) are the geometry matrix, the observation vector and the weight matrix for the subset \(p_{k}\). \(\Sigma_{p_{k}}^{0}\) is a diagonal observation weights matrix \(diag\big{(}1/\sigma_{1},1/\sigma_{2},\ldots,1/\sigma_{n_{sub}}\big{)}\) and initial weights are calculated by using following relation [34]. \[\sigma_{l}^{2}=\frac{\lambda^{2}}{2\big{(}T\cdot 10^{(CNo)/10}\big{)}^{2}}+ \frac{\lambda^{2}}{4T\cdot 10^{(CNo)/10}}\] (25) \(T\) is the predetection integration time and \(\lambda=c/f_{code}\) is the code wavelength. \(CNo\) is the carrier-to-noise ratio for corresponding satellite signal.
3. Compute the residual vector \(\mathbf{r}_{p_{k}}^{0}\), scale factor \(s_{p_{k}}^{0}\) and standardized residuals, \[\begin{array}{c}\mathbf{r}_{p_{k}}^{0}=\rho_{p_{k}}-\mathbf{H}_{p_{k}}X_{p_ {k}}^{0}\\ s_{p_{k}}^{0}=1.4826\times median\Big{(}\mathbf{r}_{p_{k}}^{0}\Big{)}\\ \widetilde{\mathbf{r}}_{s_{p_{k}}}^{0}=\frac{r_{p_{k}}^{0}}{\widetilde{r}_{p_ {k}}^{0}}\end{array}\] (26)
4. Tukey bi-weight function (Table 1) is adopted to update the weighting matrix \(\Sigma_{p_{k}}\). \[\begin{array}{c}w_{l}=1\quad\text{if}\begin{vmatrix}\widetilde{r}_{s,p_{k} }^{0}\Bigg{|}<\alpha\\ \widetilde{r}_{s,p_{k}}^{0}\Bigg{|}\geq\alpha\end{vmatrix}\] (27) Here \(\alpha\) is the tuning factor to balance the accuracy and the breakdown robustness and can be selected from the look up table in Appendix A. \(\alpha=4.658\) is selected to maintain 95% efficiency. Then updated weighting matrix will be \(\Sigma_{p_{k}}=diag\big{(}\big{\{}w_{l}:l=1,2,\ldots,n_{sub}\big{\}}\big{)}\).
5. Use Equation (24) again with updated weight matrix to compute the unknowns, then we get the estimate \(X_{p_{k}}\) for the subsample set \(p_{k}\).