GNSS Positioning Accuracy Enhancement Based on Robust Statistical MM Estimation Theory for Ground Vehicles in Challenging Environments

Global Navigation Satellite System (GNSS) is the most reliable navigation system for location-based applications where accuracy and consistency is an essential requirement. The LSE (least squares estimator) has been used since the start of GNSS for position estimation. However; LSE is affected by outliers and errors in GNSS measurements and results in wrong user position. In this paper; we proposed a novel three-phase estimator for enhancing GNSS positioning accuracy in the presence of outliers and errors; relying upon the robust MM estimation theory. In the first phase; a subsampling process is proposed on available observations. IRWLS (iterative reweighted LS) is applied to all subsamples up to a predefined number of observations to obtain a positioning estimate and a scale factor. Secondly; IRWLS is applied up to the convergence point on a set of selected subsamples. The third phase involves the selection of optimum positioning solution having minimum scale factor. An outlier detection and exclusion process is applied on a probabilistic set of outlying observations to maintain the integrity and reliability of the position. Multiple simulated and real scenarios are tested. Results show high accuracy and reliability of the proposed algorithm in challenging environments.


Introduction
Location-based applications are progressively getting more traction in our daily life, including applications such as transportation, communication, precise farming, and life-rescue operations.Positioning accuracy and consistency are the key requirements for the applications such as autonomous land vehicle driving, unmanned aerial vehicle (UAVs) and in life-critical applications.Global Navigation Satellite System (GNSS) is one of the most trustworthy positioning systems.The establishment of new GNSS constellations like the BeiDou and the Galileo has greatly increased the number of visible satellites around the global as well as new frequency signals and satellites are also added in existing GNSS constellations.These developments make it possible to use signals from multiple constellations to enhance the accuracy of user position.However, a key challenge for multiconstellation GNSS positioning is to identify the trustable satellite measurements from all visible satellites for getting true position.The quality of the estimated position depends on the ability of the receiver on getting position in the presence of outlying observations and on the receiver environments.Multipath and NLoS (non-line of sight) signals can cause errors in the satellite observations for the receiver moving in challenging environments such as urban canyons, under the elevated highway or under the dense trees.Hence it is a problematic task to achieve the positioning accuracy and reliability in such environments and it is required to adopt the positioning estimation techniques to prevent the effect of possible wrong or outlier satellite observation on the user position.
A number of algorithms and methods have been proposed to improve the GNSS receiver performance in terms of the positioning accuracy.The geometric dilution of precision-(GDOP) based best satellite selection is one of the proposed methods [1][2][3][4][5].The fundamental idea is to select a subset of observations having a minimum GDOP factor, though the rejection of a good observation can also lead to the minimized GDOP.Receiver autonomous integrity monitoring (RAIM) is another approach to reject outlier observations using pseudorange measurement redundancy to maintain the positioning accuracy by isolating the contaminated measurements from the good ones [6][7][8], and is primarily designed for aviation applications [9,10].However, for user position estimation in urban environments, visible and available satellites are rapidly changing and it is difficult to use the RAIM for satellite/pseudorange selections.Besides the satellite/observation selection techniques, several researchers also focused to improve the estimation algorithms to enhance the positioning accuracy.Least-square estimators (LSE) has been used as the fundamental algorithm for GNSS positioning since the start [11].LSE estimates the position of user by making a relationship between the unknown parameters and satellite to receiver geometry.However, it is assumed during LSE-based estimation, that the variance-covariance matrix of the observations, also called the dispersion matrix, is an identity matrix i.e., all the observations have same accuracy and are uncorrelated.On the other hand, real time GNSS observations have different level of accuracies depending upon the satellite elevation angle, signal strength and other radio interferences.Multiple algorithms and methods have been proposed to improve the GNSS positioning accuracy.WLSE (weighted least-squares estimation) schemes are proposed for calculating and assigning weights to the GNSS observations.In [12], authors proposed a WLSE (weighted least-squares estimator) where signal to noise ratio and observation variance are the weight model parameters.Tabatabaei [13] proposed a fuzzy weighted least squares method to enhance the GPS and the GLONASS (RUSSIAN GNSS) integrated navigation solution and 12 rules were defined based on the GDOP and the elevation angles of satellites for assigning weights to the GNSS observations.In [14], authors studied, tested and compared multiple weighing schemes for GNSS positioning relying upon the elevation angle and signal to noise ratio.In [15], two GNSS observation weighing models, sigma-elevation and CN0-elevation, are studied and compared by Sarab Tay.For the WLSE algorithms, satellite elevation angle and SNR (signal to noise ratio) are the primary focus of the researchers to compute the observation weights.However, relying on the elevation angle is not always useful because it is possible for the signal, from high elevation satellite, to be disrupted from the multipath and radio interferences especially in urban canyons.The concept of non-linear statistics (or robust statistics) has been also studied in literature by several researchers to improve the GNSS positioning accuracy.In [16], neural networks are used to improve the performance of least absolute and least-squares estimators for robust GNSS positioning.However, the primary disadvantage is the computational cost in addition to the outfitting problem.In [17], the concept of robust M-estimator is used to improve the performance for the GNSS positioning by calculating and adding the weight matrix during the innovation step.In [18] multiple robust statistical estimation techniques have been discussed including M-estimators, L1 estimators and genetic algorithm.The proposed estimation process in [18] determines a set of all possible solutions by using the minimum required number of satellites (e.g., in case of GPS, four satellites are required for positioning) using LSE.The median of remaining estimates is considered as the user position after removing possible wrong estimates and a genetic algorithm is used to minimize the estimation error.The fundamental problem for this scheme is the number of possible solutions.The addition of new constellations increases the number of available satellites for positioning in a single epoch and the number of solutions as well with the increase in computational cost.Secondly, the positioning process starts with an unknown random initial estimate such as a least-squares estimate and the iterative reweighting scheme can ignore perfectly good measurements resulting to a wrong convergence point with a bad initial estimate.Some researchers focused on the development of techniques to minimize the effect of noise parameters, such as ionospheric, multipath etc.In [19], authors introduced the concept of the adaptive carrier smoothing to attain an optimal carrier smoothing time to minimize the cost functions related to white noise multipath and correlated multipath errors.In [20], authors studied the application of DGNSS (differential GNSS) in smartphones for enhancement of the GNSS positioning accuracy in challenged environments.However the primary drawback for DGNSS-aided navigation is that a base station must be required for transmission of useful information and without this information host system is not able to maintain the required accuracy.
In this paper, we present a novel efficient estimation process for the GNSS positioning using signal data from multiple constellations.The proposed idea is based on the nonlinear statistical estimation theory.The proposed estimator is derived from the robust MM class estimators, proposed and named by Yohai [21].The MM estimator is a class of robust statistics for more efficient statistical data analysis with a high breakdown point.Breakdown point is the proportion measure of the outliers that can be treated efficiently without spoiling the results, and for the MM estimator, the breakdown point is 50%.For a reliable initial estimate, a weighing scheme is provided using prior knowledge of signal-to-noise ratio and satellite geometry of all visible satellites.A set of available pseudorange observations is divided into small subsets with adaptive subset size.A three-step IRWLS estimation process is applied to the individual subsample to minimize the effect of observation errors on the positioning accuracy.An outlier identification and exclusion process is proposed on the set of probabilistic outlying observations to maintain the reliability and the integrity of the positioning solution.Moreover, this paper is also presented as a concise tutorial for a better understanding of the proposed estimation algorithm in both theoretical and mathematical aspects.
The rest of the paper is organized as follows.Section 2 recalls the idea of the robust estimation theory and proposed robust estimators are reviewed.Section 3 introduces the theory and the application of the LSE model for GNSS positioning and details of the classic outlier detection exclusion process.In Section 4, mathematical and theoretical aspects of the robust MM estimation theory are presented.Details of the custom designed software receiver for performance analysis is presented in Section 5.The proposed estimation algorithm for multiconstellation positioning is explained in Section 6. Section 7 deals with the outlier detection and exclusion process in detail.In Section 8, a detailed discussion is provided related to key modifications in the original statistical estimation procedure and a very similar GNSS positioning scheme [18], based on the robust statistics, is analyzed and compared with our proposed estimation scheme.Results of the proposed estimation algorithm are analyzed and discussed against LSE and WLSE algorithms for several field and simulation tests in multiple scenarios and environments in Section 9. A conclusion is drawn in the last section with future enhancements.

Review of Robust Statistics
Robust estimation procedures have been developed and studied in statistical data analysis to replace the LSE for unknown estimation.Robust estimation algorithms perform statistical analysis of test data to detect and exclude the outlier data measurements for more accurate estimation.In normal Gaussian distribution conditions, LSE has several advantages such as simplicity of calculation.In addition, the properties of the stochastic and functional models do not change from start to end.
However, multiple classes of robust estimators have been introduced for data sets having multiple outlying observations and with the nonlinear behavior.A robust estimation algorithm is defined on three features: efficiency, stability and breakdown point.Efficiency is the performance of an estimator for an assumed model.Stability refers to the size of impairment in the performance for small deviations from the model assumptions.A robust estimation process must be stable in the sense that small deviations from the model assumptions should not impair the performance.The measure of the maximal fraction of outlier observations in any data for a robust estimator, without spoiling the final estimate, is defined as the breakdown point.A robust procedure should have an optimal or near optimal efficiency at the assumed model.It should be stable in case of small deviations from the assumed model and avoid catastrophe in case of large deviations from the model.The Least Median Square Estimator (LMSE), proposed by Hampel [22] and further improved by Rousseeuw [23], is the first estimator that achieved the maximum asymptotic break down point of 1/2.Conventional LSE is based on minimizing the sum of the squared residuals.In the LMSE, the sum is replaced by the median of the squared residuals.However, the relative efficiency of this estimator is 0 with respect to LSE because of its low rate of convergence [24].In [25], Rousseeuw and Yohai proposed the first high breakdown robust regression estimator-named S-estimator-that can achieve both high breakdown point and convergence rate.However, Hössjer proved that the regression S-estimator has maximum asymptotic efficiency 0.329 at the normal distribution [26].Under the condition of normality, Rousseeuw introduced the concept of weighted LSE (WLSE) by designing a procedure to minimize the weight of outlying observations to reduce the effect of the outliers [27].Later on, this concept became very popular in navigation research where weights are assigned based on the geometry (elevation and azimuth) and quality of signals [28,29].Yohai proposed a robust estimator that can achieve nearly optimal efficiency as well as maximum breakdown point and called as MM-estimator [21].Gervini and Yohai introduced a robust and efficient reweighted least squaress estimator (REWLSE), which combines a high breakdown point with the full asymptotic efficiency for normal errors [30].

Least Squares Estimation Model for GNSS Positioning and Outlier Correction
LSE has been used for the GNSS positioning since the start and still is the key part of many GNSS positioning algorithms such as WLSE, IRLS, etc.A standalone LSE algorithm is not a trustworthy estimator especially in GNSS-challenging environments and this is proved by experimental results in Section 9.In this section, we provide a theoretical understanding of the LSE for GNSS positioning because all positioning algorithms follow the same or a very similar approach.

Least Squares Position Estimation
In the classical GNSS-based position estimation process, LSE recursively calculates the offset of the receiver's position relative to the linearization point by using measured pseudorange errors.Let ∆ρ denotes the pseudorange error by taking the difference between the measured pseudorange ρ and the predicted pseudorange ρ: ρ is the geometric distance between the estimated receiver position and the satellite position: x s l , y s l , z s l are the position coordinates of the l_th satellite, [x u , y u , z u ] are the estimated receiver's coordinates and dt clk is the receiver clock correction.c is the speed of light.The relation between ∆ρ and ∆X u (difference between the predicted Xu and the true X u receiver positions) is expressed as, H is the geometry matrix and ξ is the measurement noise vector.
Xi u is the estimated position for i-th, Xu i−1 is the estimated position in last iteration and ∆X u i is the new estimated position error.

Receiver Autonomous Integrity Monitoring Model for Outlier Detection and Removal
RAIM is a self-monitoring process in a GNSS receiver to provide protection against the measurement faults.The classic RAIM scheme uses observation residual vector from the LSE for outlier prediction and tracking.
r ρ is the observation error residual vector.Replacing ∆X u from Equation (4) into Equation ( 6), The test statistic for fault isolation is defined as the sum of the squared residual errors, Under a fault-free positioning condition, the corresponding test statistic υ follows a chi-square (χ 2 ) distribution with n − ν degree of freedom where n is the number of measurements and ν is the number of unknowns, and should be less than a predefined threshold T th .The system false alarm probability P FA and misdetection probability P md are used to calculate the T th .
P md = P{ υ < T th |fault in sat} (10) Mathematically P FA and P md are expressed as [6], Γ(α) and I a (t) are the gamma and the first order Bessel function respectively; σ 2 is the variance of measurement error and s If the test statistic exceeds the threshold, the presence of fault is confirmed.Referring to [7,8], the standardized outlier test is defined as, e l = [0, . . . , 1,. . . 0]T , the l_th element equal to 1 while others are zeros.Q v is the cofactor matrix of the estimated residuals and is defined as, Appl.Sci.2018, 8, 876 6 of 23 A second test is performed to isolate the faulty measurements.

MM Estimation Theory for Positioning
Conventional statistical estimation schemes are based upon some assumptions about the observations.If the real observations deviate from the assumptions, estimation might deviate from the truth, such as LSE.The robust estimator tries to dynamically search a best fitting criterion to the real observations to reduce the effect of unusual data.Generally, a robust estimator minimizes the cost function J(r), ρ l is the l_th observation, X is the estimated unknown (here we consider the unknown term as generic), H l is the transformation vector for the l_th observation, and r l = ρ l − H T l X is the error residual for the l_th observation.µ(r) is the objective function that gives the contribution of each residual to the cost function J(r).For the least squares estimator µ(r) = r 2 , while in the robust estimation process, a scaled observation residual is minimized.The scaled estimation process minimizes the weight of out of the bound observations dynamically, s is the scale factor.r s is the scaled residual r/s.Let if ψ be the derivative of the objective function for the unknowns X, ψ = µ | X .Differentiating the objective function in Equation ( 17) and setting the derivative to 0, By introducing a weighting function w r s,l to consider the effects of each scaled residual r s,l , Equation ( 18) can be modified as, Solution of Equation ( 19) is a weighted least-squares problem.However, the weights depend upon the scaled residuals, the scaled residuals depend upon the estimated coefficients, and the estimated coefficients depend upon the weights.Hence, it is an iterative process, iteratively reweighted least squaress (IRWLS) [31,32].Three most popular objective functions and their corresponding weighting functions are shown in Table 1.

Name
Objective Function Weight Function r is the observation residual and it can be either true or scaled.The α, for Huber and Bi-square, is a tuning constant with optimal value to attain more resistance against outliers (high breakdown point) as well as the convergence efficiency.α is a tradeoff parameter in terms of robustness and efficiency.The table of α is provided in Appendix A to represent the relation between efficiency, robustness and tuning constant.α is selected such that to balance the performance of high breakdown point and efficiency.An outline of robust estimation procedure is defined in the following steps: 1.
Find an initial estimate, X 0 with IRWLS (Appendix B) or subsampling method (Appendix C).
In IRWLS (simplest subsampling) procedure, all observations are used in X 0 estimation.However, in subsampling procedure, subsamples (subsets) of observations are selected each with size of n sub , number of observations in one subset.For each subsample, a solution is computed using IRWLS.Best candidate solution is selected as X 0 and the initial residuals are estimated, r l (X 0 ) = ρ l − H l T X 0 for each observation.Then we get a residual vector Multiple scale estimate functions are defined, however the following median function is employed to compute the initial scale value s 0 [33],

3.
Compute initial scaled residual vector r s 0 = r 0 /s 0 .Use the weighting function w( r s ) of one method, Huber or Bi-square, shown in Table 1 to compute corresponding weights w l for each observation.

4.
Solve the estimated unknown X i at the i-th iteration (i = 1, 2, 3, . . . ) by using the weighted least-squared method, H is the transform matrix composed of H l , and Σ is the weighting matrix with the weights w l (l = 1, 2, . . ., n) as diagonal elements.

5.
Update residuals r l (X i ) = ρ l − H T l X i with the new estimate X i , calculated in last step.Then compute scaled residual r s,l = r l (X i )/s i−1 .Note that we still use last-step's scale factor to get new scaled residual because updated scale factor can only be obtained later.6.
Use the objective function µ( r s ) to compute transformed scaled residual r µ s,l = µ( r s,l ), and get a transformed scaled residual vector r Update scale factor s i with Equation ( 20) but with the transformed scaled residual vector r µ i .After that go back to step 4 to start a new iteration.
The first three steps are only used for initial estimate, and the last four steps are iteratively used for the later estimates.

GNSS Receiver Positioning System Architecture
The complete custom designed GNSS software receiver architecture is shown in Figure 1.The first step comprises to receive original carrier frequency GNSS signals using a wideband antenna and signals are down-converted to IF (intermediate frequency) data for further processing in software receiver.Receiver signal tracking module processes the IF data and performs two tasks: first the satellites, originating the signal data, are acquired/identified, secondly the satellites are tracked using the stream of signals.The next step involves the transformation of signal tracking parameters to the satellite observations and navigation information.This information is used by estimation algorithm to calculate the receiver location information i.e., position, velocity, time etc.
In Table 2, specifications of the custom designed software receiver are presented.For performance comparison, LSE, WLSE and the proposed estimator are implemented under same receiver specifications.Non-coherent scalar architecture is used for tracking satellites with the early-minus-late discriminator for code phase and cross discriminator for Doppler frequency.CNo and other observation parameters are averaged over 1 second.To reduce the dynamic stress errors, 3rd order phase lock loop (PLL) assisted by frequency lock loop (FLL) is used.First order DLL (delay lock loop) is used for code phase tracking.For n + m visible satellites from BeiDou and GPS constellations, tracking module tracks the satellite signals using Doppler frequency and code phase ( f , φ).Doppler frequency and Code phase information is used to estimate the satellite to receiver pseudorange rate .ρ and the satellite pseudorange observations (ρ) respectively.This information along other satellite information is transferred to position estimation module.Estimation process computes the receiver position in addition to maintain the integrity and reliability of calculated position.
Appl.Sci.2018, 8, x FOR PEER REVIEW 8 of 23 receiver.Receiver signal tracking module processes the IF data and performs two tasks: first the satellites, originating the signal data, are acquired/identified, secondly the satellites are tracked using the stream of signals.The next step involves the transformation of signal tracking parameters to the satellite observations and navigation information.This information is used by estimation algorithm to calculate the receiver location information i.e., position, velocity, time etc.In Table 2, specifications of the custom designed software receiver are presented.For performance comparison, LSE, WLSE and the proposed estimator are implemented under same receiver specifications.Non-coherent scalar architecture is used for tracking satellites with the earlyminus-late discriminator for code phase and cross discriminator for Doppler frequency.CNo and other observation parameters are averaged over 1 second.To reduce the dynamic stress errors, 3rd order phase lock loop (PLL) assisted by frequency lock loop (FLL) is used.First order DLL (delay lock loop) is used for code phase tracking.For nm  visible satellites from BeiDou and GPS constellations, tracking module tracks the satellite signals using Doppler frequency and code phase   , f  .Doppler frequency and Code phase information is used to estimate the satellite to receiver pseudorange rate    and the satellite pseudorange observations (  ) respectively.This information along other satellite information is transferred to position estimation module.Estimation process computes the receiver position in addition to maintain the integrity and reliability of calculated position.In addition to the estimation process, an outlier identification and exclusion procedure is also introduced for verification of outlier free positioning.A reliability test is used to check whether the outlier observations is present or not.If the presence of an outlier observation is confirmed,  In addition to the estimation process, an outlier identification and exclusion procedure is also introduced for verification of outlier free positioning.A reliability test is used to check whether the outlier observations is present or not.If the presence of an outlier observation is confirmed, probabilistic outlier observations are approximated.Selected observations are analyzed for the outlier selection and the corresponding observation is removed.

Proposed Multiconstellation Positioning Algorithm
The proposed positioning algorithm calculates the navigation information i.e., position, velocity etc. of the object where antenna is fixed.In this paper, BeiDou and GPS satellite observations are used in combination as a multiconstellation GNSS system.Let n be the number of observations from the first constellation and m be the number of observations from the second constellation (first can either be any GPS or BeiDou and other is the second).The first step involves the selection of subsets.In classic subsampling based estimation schemes, for a statistical data analysis, a small number of subsamples are drawn in random to process from all possible combinations and each subsample contains a minimum required number of observations.Subsamples are drawn in random because of large number of data samples and is not possible to process all combinations.In contrast, GNSS observations are finite and by using a random technique, it is possible that a good observation is thrown away.Hence it is more appropriate to consider all observations.If there are N number of combinations from all n + m observations with n sub observations in one combination, there will be N possible estimates G = X1 , X2 , . . ., XN .In this process, n sub is an essential parameter that needs to be chosen carefully because the position estimation time is directly related to number of subsamples.For example, if there are 20 satellite observations are available from two constellations and 5 unknowns (three positions, two clock corrections) are required to be estimated, minimum 5 observations are required for estimation of unknowns and there will be 15,504 number of combinations and it is not possible to process all combinations in the GNSS receiver.Using hit and trial approach, it is found that it is possible to relate the sample size n sub to CNo such that the n sub = (n + m) − k obs CNo<thr .k is the number of observations having CNo smaller than a defined threshold.
Because the presence of contaminated observations, the final position estimate Xrb must meet the following condition, s rb is the corresponding scale factor to the robust estimate Xrb .The proposed estimation procedure attains the required accuracy and consistency in two stages.The first stage is aimed to compute a reliable initial estimate.The second stage tries to attain the required accuracy.The two-stage algorithm is described in the following.

1.
Select P subsamples of size n sub from all n + m observations.The number of the subsamples is, 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 + , here 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 = p 1 , p 2 , . . ., p N .We compute the initial estimate X 0 p k with the WLSE method, Here H p k , ρ p k and Σ 0 p k are the geometry matrix, the observation vector and the weight matrix for the subset p k .Σ 0 p k is a diagonal observation weights matrix diag 1/σ 1 , 1/σ 2 , . . ., 1/σ n sub and initial weights are calculated by using following relation [34]. (25) T is the predetection integration time and λ = c/ f code is the code wavelength.CNo is the carrier-to-noise ratio for corresponding satellite signal.

3.
Compute the residual vector r 0 p k , scale factor s 0 p k and standardized residuals, Tukey bi-weight function (Table 1) is adopted to update the weighting matrix Σ p k .
Here α is the tuning factor to balance the accuracy and the breakdown robustness and can be selected from the look up table in Appendix A. α = 4.658 is selected to maintain 95% efficiency.Then updated weighting matrix will be Σ p k = diag w l : l = 1, 2, . . ., n sub .

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 .
Apply steps 1∼5 to all subsets and we get a set of subsample estimates and their corresponding scale factors, Next sort up the scale vector s in ascending order as s asc .Take the 5 minimal scale factors from s asc and their corresponding estimates (the number 5 is determined empirically by the authors, it can be changed).Up to now we obtained the most robust 5 initial estimates β = {X 1 , . . . ,X 5 } asc , their relative scale factors s = {s 1 , . . . ,s 5 } asc and corresponding subsamples P = p X 1 , . . .p X 5 .pX q Appl.Sci.2018, 8, 876 is the subsample used to compute X q where X q ∈ β .We apply the above steps 1∼5 iteratively to the selected subsets P. The iteration stops only when the convergence point is achieved, k cc is set as 0.001.The only difference in this step is the weighting function that we adopt Tukey bi-weight function type 2, Once all the selected 5 subset estimates have got converged, the one with minimal scale factor will be taken as the most robust initial estimate X 0 rb for the final robust estimation.

2.
Based on the robust initial estimate X 0 rb and its scale factor s o rb , we use all n + m observations at this step to get the final robust estimate.The computation procedure is the same as step 1∼5, or equivalently as the IRWLS algorithm in Appendix B. This iteration will not stop until the convergence point, Equation ( 29) is met, and the final estimate X rb is obtained.

Outlier Identification and Exclusion
The RAIM is a commonly used technique in the GNSS receivers for assurance of integrity and reliability of the estimated position.The fundamental task of the RAIM is the identification and exclusion of outlying observations for error free positioning.In the classic RAIM approach, the least-squares residual is used as the test statistic for outlier identification.As a rule of thumb, if the sum of squared residuals follows the chi-squared distribution, presence of outlier is rejected, otherwise observations are verified for outlier and the observation with the highest residual error is considered as the outlier.The primary drawback of this approach is that it is also possible that a wrong estimated position or exclusion of a wrong observation may lead to successful chi-square test.Conversely, in this paper, our proposed estimation technique itself maintains the integrity of estimated position by selection of the best fit observations.If ε proportion of the observations are outliers, the probability of the outlier free observation is µ = (1 − ε) n sub , n sub is the number of observations in one subsample and the probability of at least one good subsample is 1 − (1 − µ) N , N is the total number of subsamples.
However, to enhance the positioning accuracy and reliability, we also proposed a two-step outlier detection and exclusion process.The proposed outlier detection and exclusion scheme is explained in the following,

1.
A predetection test is performed to predict the possible outlier observations prior to the FDE (fault detection and exclusion) process.A predefined threshold level is used to predict the outlier observations based on the confidence level and on the confidence interval.A confidence interval is the range of values those act as the good estimates and the confidence level is the proportion of the confidence interval that contains the true value of their corresponding parameter.The confidence interval can be calculated from data [35,36].However, the desired confidence level is defined by the user.Let cl is the confidence level in predetection test, a threshold th cl can be computed using normal inverse cumulative distribution function (also called quantile function).
Here κ = (cl + 1)/2 is the probability for corresponding confidence level with mean µ and standard deviation σ and F t cl µ , σ = dt is the cumulative distribution function (CDF).The scaled observation residual is computed as, Here Y = ( ρ − ρ)/s sc and S = I − H H T H −1 H T .As soon as γ = n+m ∑ l=1 u 2 sc,l follows χ 2 distribution with n − ν degree of freedom, the estimated position is considered as reliable, where is the number of observations used for integrity monitoring and ν is the number of unknown to be estimated.

2.
Otherwise, a set of observations is selected as the potential outliers based on the following relation, Let ϕ = arg u sc,l >t cl ρ outlier,k 1 , ρ outlier,k 2 , . . ., ρ outlier,k n is the set of predicted outlier observations.k l is corresponding observation number.The FDE on the predicted set of outliers is applied as follows, Any faulted observation will satisfy the ε j ≥ |ε l |∀l .

Discussion
Here we study some key modifications and properties of the proposed algorithm.As described in Section 6, the fundamental algorithm is a two-stage process; the first stage involves the initialization of a best fit estimate in the presence of outliers while the second stage involves the convergence of the initial estimate.Subsampling process is the key point to find an initial robust estimate.Subsampling is also the most critical task as the number of subsamples determines the computational cost.Thus the number of observations per subsample is the most significant factor in determining the computational load.In our proposed scheme, n sub is a selectable factor, hence it is possible to tune the computational cost of algorithm.
The LSE estimator has low complexity level in terms of computational cost.IRWLS contains more complex structure than simple LSE as it requires additional arithmetic operations related to residual evaluation and updating of the observation weights and these operations are performed in iterative manner.In statistics, a huge amount of observation data, either qualitative or quantitative, is needed to be analyzed before making final decision.There are numerous statistical estimation techniques, such as robust and multivariate regression procedures [31,32], those provide general interpretation of data based on random samples drawn from original dataset as it is not possible to process all datasets, and additionally enough processing power is available with no real time requirements.On the other hand, the dataset such as GNSS has very small number of observations with certain characteristics, such as GDOP, and random sampling technique can lead to a wrong estimate.A multiconstellation system can provide better reliability and integrity, thus each subsample should contain the observations from all available constellations.And the n sub factor is managed in such a way that at least one subsample must contain all best fit satellites from both constellations.
Several authors proposed similar approaches for GPS positioning based on robust statistics.The idea in [18], is discussed in details to compare with our proposed method as there are some similarities among two algorithms.The key idea, in [18], is to find all possible solutions based on the minimum number of required observations per solution, and a genetic algorithm is applied for outlier detection and exclusion.CNo of GPS and BeiDou visible satellites is shown in Figure 2 for opensky data.similarities among two algorithms.The key idea, in [18], is to find all possible solutions based on the minimum number of required observations per solution, and a genetic algorithm is applied for outlier detection and exclusion.CNo of GPS and BeiDou visible satellites is shown in Figure 2 for opensky data.Signals from all visible satellites have well CNo .Thus it is hard for any observation to be an outlier and applying subsampling theorem in such a scenario is an impractical approach.And multi-GNSS positioning further increases the complexity in using minimum number of required number of observations.On the minimum requirement criteria basis, there will be multiple subgroups those will contain either GPS only or BeiDou only observations, thus leaving the advantages of multiple constellation positioning.

Performance Validation and Analysis
Multiple scenarios and environmental circumstances are assessed to validate the performance of the proposed estimator.Considered scenarios comprise static as well as dynamic environments.Two static locations are rooftop (for line of sight signal propagation model) and a location surrounded by skyscrapers (for multipath in urban environment) in the downtown of shanghai.In Figure 3a  Signals from all visible satellites have well CNo.Thus it is hard for any observation to be an outlier and applying subsampling theorem in such a scenario is an impractical approach.And multi-GNSS positioning further increases the complexity in using minimum number of required number of observations.On the minimum requirement criteria basis, there will be multiple subgroups those will contain either GPS only or BeiDou only observations, thus leaving the advantages of multiple constellation positioning.

Performance Validation and Analysis
Multiple scenarios and environmental circumstances are assessed to validate the performance of the proposed estimator.Considered scenarios comprise static as well as dynamic environments.Two static locations are rooftop (for line of sight signal propagation model) and a location surrounded by skyscrapers (for multipath in urban environment) in the downtown of shanghai.In Figure 3a,b similarities among two algorithms.The key idea, in [18], is to find all possible solutions based on the minimum number of required observations per solution, and a genetic algorithm is applied for outlier detection and exclusion.CNo of GPS and BeiDou visible satellites is shown in Figure 2 for opensky data.Signals from all visible satellites have well CNo .Thus it is hard for any observation to be an outlier and applying subsampling theorem in such a scenario is an impractical approach.And multi-GNSS positioning further increases the complexity in using minimum number of required number of observations.On the minimum requirement criteria basis, there will be multiple subgroups those will contain either GPS only or BeiDou only observations, thus leaving the advantages of multiple constellation positioning.

Performance Validation and Analysis
Multiple scenarios and environmental circumstances are assessed to validate the performance of the proposed estimator.Considered scenarios comprise static as well as dynamic environments.Two static locations are rooftop (for line of sight signal propagation model) and a location surrounded by skyscrapers (for multipath in urban environment) in the downtown of shanghai.In Figure 3a  For stationary scenarios, accuracy is examined in terms of 90% position error radius (RCEP = radial circular error position) and position error mean distance (PMD) from the center.Results are analyzed against LSE and WLSE (weighted least squares estimation) algorithms.In this paper, we applied the SIGMA − ε weight model [37] for observation weight estimation.
σ l (i) is the variance of lth observation at ith epoch.a is the scaling constant and depends upon the noise characteristics of the tracking loop or carrier tracking loop noise bandwidth.And the WLSE is defined as; Σ −1 is the weight matrix where Σ is the diagonal observation variance-covariance matrix.

Open Sky Scenario
In the open sky scenario, the only errors are Gaussian white noise from tracking loops, atmospheric and receiver/satellite clock corrections.In Table 3, GNSS signal information for open sky scenario is represented.
Horizontal radial error for open sky environment is shown in the Figure 4 in east and north error terms.There is a very little difference between RCEP of three models.However, PMD of the proposed scheme is larger than the other two schemes (except the difference is substantively small).The primary cause of better performance of LSE/WLSE is the presence of Gaussian white noise.As RCEP provides the 90% position error radius, results indicate that WLSE and the proposed schemes are more robust although the mean position error distance of LSE is smaller.Two additional simulations are also created by modifying the noise characteristics of observations.The statistics of the first emulated scenario are presented in the following Table 4.This scenario is designed to evaluate the positioning accuracy in the presence of fixed outliers.Measurement errors stay fixed for a satellite over a defined period of time.The measurement error is kept small because a large error is easy to detect by a simple algorithm however it is more difficult to detect if the error is not very large and show significant effect.Results, in Figure 5, showed that the performance of our  The statistics of the first emulated scenario are presented in the following Table 4.This scenario is designed to evaluate the positioning accuracy in the presence of fixed outliers.Measurement errors stay fixed for a satellite over a defined period of time.The measurement error is kept small because a large error is easy to detect by a simple algorithm however it is more difficult to detect if the error is not very large and show significant effect.Results, in Figure 5, showed that the performance of our proposed estimation scheme is better than LSE and WLSE in both RCEP and PMD.The statistics of the first emulated scenario are presented in the following Table 4.This scenario is designed to evaluate the positioning accuracy in the presence of fixed outliers.Measurement errors stay fixed for a satellite over a defined period of time.The measurement error is kept small because a large error is easy to detect by a simple algorithm however it is more difficult to detect if the error is not very large and show significant effect.Results, in Figure 5, showed that the performance of our proposed estimation scheme is better than LSE and WLSE in both RCEP and PMD.

Emulated Scenario 2
The objective of this scenario is to evaluate the positioning reliability and accuracy for unstable observation data.Figure 6 represents the dynamic environment simulation scenario.

Emulated Scenario 2
The objective of this scenario is to evaluate the positioning reliability and accuracy for unstable observation data.Figure 6 represents the dynamic environment simulation scenario.Specifications include the details about PRN, measurement error, number of outliers etc. Navigation module receives the signal data from tracking loops and errors are introduced in the real measurements dynamically based on the error specifications.PRN/Sat IDs are selected at random and are changed every 10 s. Results are presented in Figure 7. Position mean distance error of our proposed estimator is significantly smaller as well as RCEP.On the other hand, WLSE is also performed well against LSE.One important point to consider is that if the observation is effected by un-modeled noise sources such as radio interferences, the performance of WLSE is degraded though observations with high CNo because the weight estimation scheme for WLSE is based on the observation CNo .

Multipath Scenario in Urban Canyons
This scenario is designed to analyze the performance of the proposed algorithm in severe urban environments.Multipath and line of sight problem are the primary source of errors in urban Specifications include the details about PRN, measurement error, number of outliers etc. Navigation module receives the signal data from tracking loops and errors are introduced in the real measurements dynamically based on the error specifications.PRN/Sat IDs are selected at random and are changed every 10 s. Results are presented in Figure 7. Position mean distance error of our proposed estimator is significantly smaller as well as RCEP.On the other hand, WLSE is also performed well against LSE.One important point to consider is that if the observation is effected by un-modeled noise sources such as radio interferences, the performance of WLSE is degraded though observations with high CNo because the weight estimation scheme for WLSE is based on the observation CNo.

Emulated Scenario 2
The objective of this scenario is to evaluate the positioning reliability and accuracy for unstable observation data.Figure 6 represents the dynamic environment simulation scenario.Specifications include the details about PRN, measurement error, number of outliers etc. Navigation module receives the signal data from tracking loops and errors are introduced in the real measurements dynamically based on the error specifications.PRN/Sat IDs are selected at random and are changed every 10 s. Results are presented in Figure 7. Position mean distance error of our proposed estimator is significantly smaller as well as RCEP.On the other hand, WLSE is also performed well against LSE.One important point to consider is that if the observation is effected by un-modeled noise sources such as radio interferences, the performance of WLSE is degraded though observations with high CNo because the weight estimation scheme for WLSE is based on the observation CNo .

Multipath Scenario in Urban Canyons
This scenario is designed to analyze the performance of the proposed algorithm in severe urban environments.Multipath and line of sight problem are the primary source of errors in urban

Multipath Scenario in Urban Canyons
This scenario is designed to analyze the performance of the proposed algorithm in severe urban environments.Multipath and line of sight problem are the primary source of errors in urban environments.Experiment location is surrounded by high skyscrapers as shown in the Figure 1b.In Figure 8 CNo of visible satellite is shown from GPS and BeiDou satellites.
environments.Experiment location is surrounded by high skyscrapers as shown in the Figure 1b.In Figure 8 CNo of visible satellite is shown from GPS and BeiDou satellites.

High Dynamic Moving Scenario
The most critical requirement for every GNSS receiver is the performance under high dynamics when moving in the severe Examples of harsh surroundings are moving under trees, under elevated roads and in urban canyons with partial visible sky.High dynamic characterizes the  Positioning results are provided in the Figure 9 in terms of RCEP and PMD.As RCEP represents the accuracy of estimated position during most of time, it can be found in the figure that the proposed estimation technique significantly improved the positioning accuracy.
environments.Experiment location is surrounded by high skyscrapers as shown in the Figure 1b.In Figure 8 CNo of visible satellite is shown from GPS and BeiDou satellites.

High Dynamic Moving Scenario
The most critical requirement for every GNSS receiver is the performance under high dynamics when moving in the severe surroundings.Examples of harsh surroundings are moving under trees, under elevated roads and in urban canyons with partial visible sky.High dynamic characterizes the

High Dynamic Moving Scenario
The most critical requirement for every GNSS receiver is the performance under high dynamics when moving in the severe surroundings.Examples of harsh surroundings are moving under trees, under elevated roads and in urban canyons with partial visible sky.High dynamic characterizes the condition when the Doppler shift may vary in the range of ±60 kHz.Although this is difficult to achieve in normal urban canyons, still at carrier-to-noise ratio of the order of 30-40 dB-Hz and sudden loss of tracking satellites make it challenging to maintain the navigation accuracy.The following scenario is designed to evaluate the dynamic performance of the proposed estimation procedure against classic LSE and WLSE algorithms.GNSS signal data is recorded in remote rural area near shanghai during driving on the road with antenna fixed on the roof of van.Ephemeris information and other navigation data is produced by offline processing of signal data.
Experiment trajectory of vehicle is shown in Figure 10.The experiment route is designed so that there must be multiple full and partial sky visibility regions.Figure 11 illustrates the CNo of GPS satellites.As the CNo of satellites represent the signal characteristics, observations from the satellites, with unstable CNo, can lead the estimator to wrong approximation.
Appl.Sci.2018, 8, x FOR PEER REVIEW 18 of 23 condition when the Doppler shift may vary in the range of 60 kHz .Although this is difficult to achieve in normal urban canyons, still at carrier-to-noise ratio of the order of 30-40 dB-Hz and sudden loss of tracking satellites make it challenging to maintain the navigation accuracy.The following scenario is designed to evaluate the dynamic performance of the proposed estimation procedure against classic LSE and WLSE algorithms.GNSS signal data is recorded in remote rural area near shanghai during driving on the road with antenna fixed on the roof of van.Ephemeris information and other navigation data is produced by offline processing of signal data.Experiment trajectory of vehicle is shown in Figure 10.The experiment route is designed so that there must be multiple full and partial sky visibility regions.Figure 11     Appl.Sci.2018, 8, x FOR PEER REVIEW 18 of 23 condition when the Doppler shift may vary in the range of 60 kHz .Although this is difficult to achieve in normal urban canyons, still at carrier-to-noise ratio of the order of 30-40 dB-Hz and sudden loss of tracking satellites make it challenging to maintain the navigation accuracy.The following scenario is designed to evaluate the dynamic performance of the proposed estimation procedure against classic LSE and WLSE algorithms.GNSS signal data is recorded in remote rural area near shanghai during driving on the road with antenna fixed on the roof of van.Ephemeris information and other navigation data is produced by offline processing of signal data.Experiment trajectory of vehicle is shown in Figure 10.The experiment route is designed so that there must be multiple full and partial sky visibility regions.Figure 11 illustrates the CNo of GPS satellites.As the CNo of satellites represent the signal characteristics, observations from the satellites, with unstable CNo , can lead the estimator to wrong approximation.However, the circumstances are quite different for NLoS (non-line of sight) conditions i.e., driving under trees.During NLoS conditions, other than Gaussian errors are dominant such as multipath, signal diffraction or reflection etc.We cannot say that the standalone positioning accuracy of the proposed scheme is very high (Figure 13b), however there is a significant improvement in positioning accuracy in comparison to LSE and WLSE.
Positioning error for moving scenario from three algorithms are analyzed in Table 5. Maximum and RMS (root mean square) position errors from the proposed estimation algorithm are significantly lower than LSE and WLSE.The maximum error are resulted when vehicle was moving under dense trees and majority of the received observations were corrupted.
Table 5. Analysis of position errors from three algorithm for moving scenario.However, the circumstances are quite different for NLoS (non-line of sight) conditions i.e., driving under trees.During NLoS conditions, other than Gaussian errors are dominant such as multipath, signal diffraction or reflection etc.We cannot say that the standalone positioning accuracy of the proposed scheme is very high (Figure 13b), however there is a significant improvement in positioning accuracy in comparison to LSE and WLSE.
Positioning error for moving scenario from three algorithms are analyzed in Table 5. Maximum and RMS (root mean square) position errors from the proposed estimation algorithm are significantly lower than LSE and WLSE.The maximum error are resulted when vehicle was moving under dense trees and majority of the received observations were corrupted.
Table 5. Analysis of position errors from three algorithm for moving scenario.However, the circumstances are quite different for NLoS (non-line of sight) conditions i.e., driving under trees.During NLoS conditions, other than Gaussian errors are dominant such as multipath, signal diffraction or reflection etc.We cannot say that the standalone positioning accuracy of the proposed scheme is very high (Figure 13b), however there is a significant improvement in positioning accuracy in comparison to LSE and WLSE.
Positioning error for moving scenario from three algorithms are analyzed in Table 5. Maximum and RMS (root mean square) position errors from the proposed estimation algorithm are significantly lower than LSE and WLSE.The maximum error are resulted when vehicle was moving under dense trees and majority of the received observations were corrupted.

Conclusions and Future Perspectives
A novel, real-time GNSS positioning estimator is presented in this paper for land vehicles in challenging environments.The proposed estimation procedure is based on robust MM estimation theory.The presence of outlier observations in GNSS measurements significantly affect the positioning results.This method is capable of approximating the user position with accuracy and reliability in the presence of outlying observations by using the most trusted observations for position estimation and excluding outliers.The proposed algorithm is compared against LSE and WLSE algorithms in terms of 90% positioning accuracy.Extensive field and simulation tests are carried out to analyze the performance of the proposed scheme in diverse environments.Results demonstrate that the proposed scheme could be a promising scheme for GNSS-based positioning with high accuracy under diverse conditions, especially in harsh and challenging environments, because LSE and WLSE are also capable to provide similar accuracy in normal outlying free environments, however they are not able to maintain the accuracy in challenging environments.For future extensions, we plan to explore the applications of the algorithm in more complex GNSS-based positioning architectures such as vector-based GNSS tracking receivers and fusion of GNSS measurements with other sensors such as IMU (Inertial Measurement Unit) and vision sensors to further improve the positioning accuracy.We also want to explore the application of the proposed estimation scheme (especially therobust initialization stage) for integration with the EKF (Extended Kalman Filter) in order to enhance the filter performance.

Figure 1 .
Figure 1.Implementation flowchart of proposed estimation scheme.

Figure 1 .
Figure 1.Implementation flowchart of proposed estimation scheme.

Figure 2 .
Figure 2.CNo of BeiDou and GPS-visible satellites for open sky data.

Figure 2 .
Figure 2. CNo of BeiDou and GPS-visible satellites for open sky data.
, the two locations are shown.In open sky scenario (roof top) all signals are line of sight with no multipath or other interferences.In contrast, there are multipath as well as other radio interferences in urban canyon scenario.Appl.Sci.2018, 8, x FOR PEER REVIEW 13 of 23

Figure 2 .
Figure 2. CNo of BeiDou and GPS-visible satellites for open sky data.

Figure 4 .
Figure 4. Horizontal radial error of open sky static scenario.

Figure 4 .
Figure 4. Horizontal radial error of open sky static scenario.

Figure 4 .
Figure 4. Horizontal radial error of open sky static scenario.

Figure 6 .
Figure 6.Flowchart of observation error simulation in true navigation data.

Figure 6 .
Figure 6.Flowchart of observation error simulation in true navigation data.

Figure 6 .
Figure 6.Flowchart of observation error simulation in true navigation data.

Figure 8
Figure 8 demonstrations that CNo of both BeiDou and GPS is unstable.And the disruption in line represents the loss of satellite visibility.Moreover the multipath signal interferences introduce additional errors in observations.Positioning results are provided in the Figure9in terms of RCEP and PMD.As RCEP represents the accuracy of estimated position during most of time, it can be found in the figure that the proposed estimation technique significantly improved the positioning accuracy.

Figure 9 .
Figure 9. HRE of proposed estimator against LSE and WLSE.

Figure 8
Figure 8 demonstrations that CNo of both BeiDou and GPS is unstable.And the disruption in line represents the loss of satellite visibility.Moreover the multipath signal interferences introduce additional errors in observations.Positioning results are provided in the Figure9in terms of RCEP and PMD.As RCEP represents the accuracy of estimated position during most of time, it can be found in the figure that the proposed estimation technique significantly improved the positioning accuracy.

Figure 8
Figure 8 demonstrations that CNo of both BeiDou and GPS is unstable.And the disruption in line represents the loss of satellite visibility.Moreover the multipath signal interferences introduce additional errors in observations.Positioning results are provided in the Figure9in terms of RCEP and PMD.As RCEP represents the accuracy of estimated position during most of time, it can be found in the figure that the proposed estimation technique significantly improved the positioning accuracy.

Figure 9 .
Figure 9. HRE of proposed estimator against LSE and WLSE.

Figure 9 .
Figure 9. HRE of proposed estimator against LSE and WLSE.
illustrates the CNo of GPS satellites.As the CNo of satellites represent the signal characteristics, observations from the satellites, with unstable CNo , can lead the estimator to wrong approximation.

Figure 10 .
Figure 10.Vehicle route for data logging (courtesy of Google Earth).

Figure 10 .
Figure 10.Vehicle route for data logging (courtesy of Google Earth).

Figure 10 .
Figure 10.Vehicle route for data logging (courtesy of Google Earth).

Figure 12a shows the 23 FigureFigure 12 .Figure 13 .
Figure 12a shows the errors in the estimated position from LSE, WLSE and our proposed algorithms.As described before, there are several open sky and partial visible sky regions, mentioned with arrows.In Figure 12b, position errors for open sky are zoomed.Estimated positions errors, for the duration of driving under trees, are showed in Figure 13.From Figure 11b, it is found that the performance of three algorithms is similar under (LoS) line-of-sight signal propagation conditions such as open sky.In LoS conditions, signals are biased with Gaussian-only white noise, hence LSE showed high positioning accuracy.Appl.Sci.2018, 8, x FOR PEER REVIEW 19 of 23

Figure 13 .
Figure 13.Estimated position errors when during driving under trees (zoomed in).

Figure 13 .
Figure 13.Estimated position errors when during driving under trees (zoomed in).
The noise might include uncompensated atmospheric delay, multipath error, satellite clock and ephemeris errors, and the Gaussian white noise.The receiver position is recursively computed,

Table 3 .
GNSS signal information for open sky scenario.

Table 5 .
Analysis of position errors from three algorithm for moving scenario.