# A new Robust Kalman filtering algorithm of unequal precision observations based on residual vectors in static precise point positioning.

ABSTRACTPrecise Point Positioning (PPP), which can achieve high-precision positioning with only a single Global Navigation Satellite System (GNSS) receiver, is a popular but challenging research topic. The traditional Robust Kalman Filter (RKF) based on innovation vectors can effectively resist outliers in certain cases. However, this filter is less effective in treating unequal precision observations for PPP, such as crowded outliers or small outliers in high-precision carrier phase observations. In this study, the residual vector is used to construct a robust factor that is sensitive to outliers. Our strategy is to apply decorrelation to this vector. Firstly, the squared Mahalanobis distances of carrier phase and pseudorange observations are used to evaluate whether measurements contain outliers in the current epoch. Secondly, the residual vector is decorrelated with respect to the residual covariance. Finally, through iteration, the robust factor for the residual vector and the gain matrix are determined, which theoretically eliminates the residual vector correlation for different observations. Our proposed modification of the RKF method has been tested using data from International GNSS Service (IGS) Stations. Results show that our RKF based on residual vectors can effectively reduce the effects of a single outlier. Given many outliers, computations must be iterated to reduce the residual vector correlation, especially for cases where more outliers can effectively be suppressed by filter divergence.

1. INTRODUCTION

PPP is a popular but difficult topic in GNSS research. Such positioning uses only a single GNSS receiver, precise orbit and clock products to achieve centimeter-decimeter, and even millimeter-level positioning (Zumberge, 1997). Both carrier phase and pseudorange measurement values are used. To reduce the influence of ionospheric delay, such ionosphere-free combination observations are used by most scholars (Kouba and Heroux, 2001; Hefty and Gerhatova, 2012; Zhang et al., 2012, 2014). In addition, the University of Calgary (UofC) PPP model has been shown by Gao and Shen (2002) to eliminate observation noise. The ionosphere-free combination of carrier phase and pseudorange observations is selected for study here, because their ionospheric delays have the same carrier frequency and they are similar in size, although they comprise positive and negative numbers, respectively. Three main ambiguity fixing methods have been widely used to overcome slow convergence and low precision in traditional PPP (Ge et al., 2008; Laurichesse et al., 2009; Collins et al., 2010; Shi and Gao, 2014). Such methods can significantly improve positioning accuracy and shorten the convergence time. Many scholars have proposed real-time PPP or PPP-RTK technology based on ambiguity fixing methods, precise orbit and clock products ( Zhang et al., 2014; Li et al., 2014, 2015a and 2015b; Dawidowicz and Krzan, 2014; Krzan and Przestrzelski, 2016).

The Kalman filter has been used extensively to estimate the state space in positioning and navigation, because the status parameters can be estimated in real time (Kalman, 1960). This filter method requires only the prior value of the status parameters, covariance, and new observation data. Yang (1997, 2001) proposed a Robust Kalman Filtering (RKF) algorithm based on innovation vectors, and an equivalent weight matrix structure with outliers. An innovation-based adaptive RKF was developed for dynamic PPP (Guo and Zhang, 2014), which introduces an equivalent covariance matrix to resist outliers, and an adaptive factor to balance the contribution of observational versus predicted information from the dynamic model system. Nie (2010) merely used RKF to resist carrier phase outliers, because their influence on positioning is stronger than pseudorange observations.

Here, we present an RKF method based on residual vectors for static PPP. To assess its effectiveness, experiments were conducted with a single outlier and multiple outliers. We also compared the results computed with traditional Kalman filter and our proposed RKF.

2. KALMAN FILTERING MODEL FOR STATIC PPP

Traditional PPP relies on both the carrier phase and the pseudorange of two unequal precision measurement values. At present, most processing strategies adopt the first-order ionosphere-free elimination of the carrier phase and pseudorange observations for such positioning, because the ionospheric delay model is more difficult to employ for estimation and introduces complexity. The combination of ionosphere-free carrier phase and pseudorange observations can be represented as (Hu et al., 2014; Xu et al., 2012):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2)

where [[phi].sub.IF] and [P.sub.IF] denote the ionosphere-free combination of carrier phase and pseudorange observations in meters, respectively. [[florin].sub.1] and [[florin].sub.2] are the frequencies of L1 and L2 GPS carrier phase signals, where [[florin].sub.1]=1 575.42 MHz and [[florin].sub.2]=1 227.60 MHz. [[phi].sub.i] and [P.sub.i] ( i = 1,2 ) are the carrier phase and pseudorange measurements, respectively. [rho] is the geometrical distance from the GPS receiver to the satellite in meters. Constant c is the speed of light in vacuum. Parameter dt represents the clock errors in seconds; the errors in the receiver and the satellite clock errors are collectively considered. Parameter T is the residual tropospheric delay after correction using the Saatamoinen model. [epsilon] and e are the noise in the ionosphere-free combination of carrier phase and pseudorange observations, respectively; the carrier phase and pseudorange observation noise occur within 0.05 and 0.5 m, respectively. Parameter [[lambde].sub.IF] is the coefficient of the ionosphere-free ambiguity solution, where [[lambde].sub.IF]=10.7 cm. Parameter N is the ionosphere-free ambiguity of the carrier phase observations. In general, satellite orbit and clock errors have decreased considerably with the use of the interpolated IGS precise orbit and clock products. In addition, the pseudorange biases are corrected using the Differential Code Biases (DCB) product from IGS. Other error source have been corrected and reduced by appropriate model or IGS products, and are not listed in the Equation (1) and (2) for simplification, such as those caused by the Earth's rotation, the satellite antenna phase center, the receiver antenna phase center, relativistic effects, tidal error, and other factors. Hardware delay error was taken into account as part of clock error and ignored the influence.

Given that the static PPP does not reflect disturbance, the influence of observation outliers alone is considered here; hence, adaptive filtering effects on positioning were not considered. In Equations (1) and (2), the estimated parameter vectors include the station coordinates, residual tropospheric delay, receiver clock error and each satellite's real ionosphere-free ambiguity parameters. The estimated parameters are expressed as:

X = [([[delta].sub.x],[[delta].sub.y],[[delta].sub.z],T, dt, [N.sup.1],... ,[N.sup.j].sup.T]

where X represents the estimated parameter vector; [[delta].sub.x],[[delta].sub.y],[[delta].sub.z] are the 3D coordinate vectors of the receiver in the ITRF08 coordinate system; j is the number of visible satellites; and T denotes the matrix transposition.

The carrier phase and pseudorange observations of the k th epoch given by the Taylor series expansion can be expressed as:

[Z.sub.k]=[H.sub.k][Z.sub.k]+[[epsilon].sub.k] (3)

where [Z.sub.k] is the observation vector that includes both carrier phase and pseudorange observations; k represents the epoch location; [H.sub.k] is the coefficient matrix of the observation; and [[epsilon].sub.k] is the observation noise vector.

The state and dynamic system equations are written as:

[X.sub.k]=[[PHI].sub.k,k-1][X.sub.k-1]+[[omega].sub.k] (4)

where [[PHI].sub.k,k-1] is the state transition matrix; here, [[PHI].sub.k,k-1] is a unit matrix because only static PPP is considered in this study, [[omega].sub.k] is the state noise vector.

Thus, the standard Kalman filtering process involves:

Forecast processing:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Update processing: ,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where X is the predicted state parameter vector; [PHI] is the state transition matrix; X is the estimated state parameter vector; [R.sub.[bar.X]] is the predicted state covariance matrix; [R.sub.X] is the posteriori covariance matrix of the state parameters; Q is the covariance matrix of the processing noise; K is the gain matrix of the Kalman filter; H is the coefficient matrix of the observation equation; P is the observation weight matrix; R is the observation noise variance matrix, such that R = [P.sup.-1] ; and V is the predicted parameter vector, also called the innovation vector.

3. RKF MODEL OF UNEQUAL PRECISION OBSERVATIONS

Two methods are often used to eliminate or reduce the impact of outliers; they involve either outliers' detection using a test statistic or robust estimation (Trasak and Stroner, 2014). The standard Kalman filter can be used only for random errors that follow a normal distribution model. If we still apply such filtering, and ignore the effect of outliers on observations, then the resultant estimated parameter will be incorrect and may lead to filter divergence. An RKF model is highly useful in this regard. At present, the RKF model includes a model for both observational and dynamic anomalies (Yang, 1997 and 2001). Here, we evaluate a new RKF model for original observation data with outliers, because the static PPP problem is considered in this study.

3.1. RKF OF ORIGINAL OBSERVATION DATA WITH OUTLIERS

The residual and innovation vectors of a standard Kalman filter can be respectively expressed as:

[V.sub.k]=[H.sub.k][X.sub.k]-[Z.sub.k] (5)

[V.sub.k]=[H.sub.k][[bar.X].sub.k]-[Z.sub.k] (6)

where [V.sub.k] is the posteriori residual vector of observation, which reflects posteriori residuals; [V.sub.k] is the predicted state residual vector of an observation, which reflects priori residuals; [[bar.X].sub.k] is the predicted state vector; [X.sub.k] is a parameter vector that is estimated based on the original observations; and [Z.sub.k] is the observation vector.

We assume that the k-1 th epoch observation vector [Z.sub.k] has no outliers, and that the k th epoch state parameter vector [[bar.X].sub.k] has a Gaussian distribution. Conversely, if the k th epoch of the original observation data contains outliers, then the observation vector has a Gaussian distribution of noise. In this case, the unknown observation parameter vector should be estimated robustly, whereas the state parameter vector is estimated using a least squares method. The extreme conditions are as follows (Yang, 2001):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)

where [P.sub.k] is the robust weight matrix of the observation vector given by:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where el denotes the elevation angle of the satellites.

Upon substituting the extreme values from Equations (3)-(6) into Equation (7), the state vector and the gain matrix can be estimated using

[[bar.X].sub.k]=[[bar.X].sub.k]+[K.sub.k]([H.sub.k][[bar.X].sub.k]-[Z.sub.k]) (8)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (9)

As per Equation (6), [V.sub.k] is affected by the observation vector [Z.sub.k], which may contain outliers. To reduce the influence of any outliers, we must construct a reasonable gain matrix [K.sub.k]. This matrix is mainly affected by [P.sub.k] and [R.sub.[bar.X]]. Because [R.sub.[bar.X]] is solved by processing noise covariance and posteriori parameters obtained from the previous epoch solution, which is unrelated to the current epoch of information. Meanwhile, the equivalent weight matrix [P.sub.k] is determined by satellite altitude angle or noise. A gain matrix can be calculated from the equivalent weight matrix. Some studies construct a robust factor directly by the gain matrix (Wang et al., 2008); thus, the size and location of any outliers is difficult to evaluate, and the gain matrix has a more complex structure. Instead, we construct the corresponding weight matrix, so that the gain matrix is not directly based on the robust factor.

3.2. CALCULATION OF THE ROBUST FACTOR BASED ON RESIDUAL VECTORS

The innovation vector is often used to construct a robust factor in traditional RKF. Here, we analyze the sensitivity of this vector and of the residual vector to outliers.

The innovation vector from Equation (6) can be expanded for carrier phase and pseudorange observations, as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (10)

where C is the non-ambiguity component of the parameter and N is the ambiguity component of the parameter.

Moreover,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [V.sub.pha] and [V.sub.pse] are the innovation vectors of the carrier phase and pseudorange observations, respectively; [H.sub.pha] and [V.sub.pse] are the coefficient matrices of the carrier phase and pseudorange observations, respectively; [Z.sub.pha] and [V.sub.pse] are the carrier phase and pseudorange observations, respectively; [H.sub.unamb] is the coefficient matrix of the non-ambiguity component; and n denotes the number of visible satellites. The carrier phase and pseudorange observations of unequal precision have the same [H.sub.unamb] matrix. In addition, we use eye(n) to denote an n-dimensional unit matrix, and zero( n) to denote an n-dimensional zero matrix.

As shown in Equation (10), coefficient matrices [H.sub.unamb] of the carrier phase and the pseudorange observations are equal, with the exception of the ambiguity coefficient; only the observation values and the ambiguity coefficient matrix differ. The value of the innovation vector can be determined using the predicted state parameter vector of the previous and current epoch observations. If the estimated state parameters in the previous epoch are accurate, then the innovation vector responds directly to outliers in the current epoch. The observational accuracy of the carrier phase is much higher than the pseudorange observations; hence, the innovation values of all satellite carrier phase vectors are approximately equal. Thus, the innovation vector can effectively identify carrier phase and pseudorange observations in the current epoch. Moreover, pseudorange observations can be used to distinguish the size and location of outliers as well. However, the innovation value of the carrier phase containing outliers will deviate significantly from the rest of the observations from other satellites, because of their more accurate in nature. Although an abnormal value is easy to detect manually, such values cannot be distinguished automatically from correct observations, because the noise error is often greater than the effect of outliers on the innovation vector. Many outliers cannot be recognized either manually or automatically because abnormal outliers severely pollute the quality of normal observation data, making these outliers indistinct. As a result, innovation-based traditional Kalman filtering is ineffective on carrier phase observations with outliers, particularly when small outlier values are present.

The residual vector from Equations (5) and (8) is estimated using a posteriori observation value; this value is used to solve the posteriori state parameters from unequal precise observations, which is determined using reasonable weights. Hence, this residual vector is more sensitive than the innovation vector to outliers in unequal precision observation data. The relationship between the residual and the innovation vectors, as well as the covariance matrix, can be expressed as (Yang, 2001):

[V.sub.k]=(I-[H.sub.k][K.sub.k])[V.sub.k] (11)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (12)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the covariance matrix of [V.sub.k]; and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the covariance matrix of [V.sub.k] . The innovation vector can be derived from the state vector, which has not been corrected by the observation vector. Therefore, the innovation vector reflects dynamic disturbance. In contrast, the residual vector is solved using the actual parameter vector, modified by the observation vector. Equation (12) indicates that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] must be less than [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] ; therefore, the residual vector is more sensitive than the innovation vector to outliers.

In summary, the innovation vector can reflect the location and size of outliers, given highly accurate observation values. However, correct observation values cannot be distinguished from values containing abnormal outliers, as the number of outliers increases. The residual vector is a posteriori state parameter that is estimated using reasonable observation weights in the current epoch; it effectively reflects both the size and location of outliers. Finally, the covariance estimation must have an adequate redundant observation value. Because this research does not have enough redundant observations to ensure statistical quality, we have constructed a robust factor based on residual vectors for comparative purposes only.

Since the residual vector is subject to a Gaussian distribution, we need to construct a density function that is similar to this distribution. The common weight functions of the robust factor are the Huber, IGG (Institute of Geodesy and Geophysics) II, and IGG III (Huber, 1981; Yang, 2002). We selected the IGG III as the Gaussian distribution function to construct the equivalent weight matrix, as follows (Yang, 2002):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (13)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] represents the equivalent weight matrix of the (i,i) element, and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is the weight matrix of the original observation noise i . [y.sub.ii] is given by:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (14)

where [k.sub.0] and [k.sub.1] correspond to a sub-parameter(ranging from 1.5 to 2 ) and an elimination value(ranging from 3 to 8 ), respectively. [s.sub.ii] is the robust factor.

Moreover (Xu et al., 2012),

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15) where [V.sub.i] represents the i th element of the residual vector, and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] denotes the covariance matrix of the residual vector.

When RKF is used on unequal precision measurements, the unequal precision of the original observation noise may result in a mismatch between the two observations in terms of [k.sub.0] and [k.sub.1]. The residuals of the carrier phase and the pseudorange observations are not completely consistent with the residual; thus, the corresponding robust factors [k.sub.0] and [k.sub.1] are discretely defined. In this study, we considered the carrier phase and the pseudorange observations to be uncorrelated; hence, the IGG III weight function was used to generate the robust factor.

Given that it is cumbersome to determine the residuals of every observation, we modified the traditional method. The Mahalanobis distance of the carrier phase and pseudorange observations for every epoch was calculated. The Mahalanobis distances reflect overall abnormalities (Chang, 2014a, 2014b and 2015) and can be used to determine whether there are outliers in the current epoch by comparing the Mahalanobis distance with a threshold value. Finally, our proposed RKF algorithm was employed. Using this procedure, processing efficiency was greatly enhanced.

The squared Mahalanobis distance of the carrier phase and pseudorange observations is expressed as:

[M.sup.2] =[V.sup.T] PV~ [[chi].sup.2]

where [M.sup.2] represents the squared Mahalanobis distance, which follows a chi-square distribution [[chi].sup.2] . According to the hypothesis testing theory, the probability of a squared Mahalanobis distance less than the quantile of the distribution [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] can be expressed as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where P(*) represents the probability. The probability that [M.sup.2] is larger than [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] should be very small, i.e., less than [alpha]. Because this strategy greatly improved the efficiency of data processing, the squared Mahalanobis distance and the quantile of the carrier phase and pseudorange observations were constructed for all unequal precision measurements, respectively.

3.3. DECORRELATION FOR RESIDUAL VECTORS

When Equation (5) is transformed, the residual vector can be expressed as:

[V.sub.k]=[H.sub.k]([[bar.X].sub.k]+[K.sub.k]([H.sub.k][[bar.X].sub.k]-[Z.sub.k]))-[Z.sub.k]

Here, the residual vector reflects the size and location of small errors. Nonetheless, the residual vectors of other satellites are affected, whenever the observations of any one satellite contain outliers. This becomes relevant, when many satellite datasets contain outliers. In such cases, constructing a robust factor and an equivalent weight matrix directly from such data is unreasonable.

Given less accurate pseudorange observations, the effect on posterior values X is negligible. Therefore, the residual vector of the carrier phase measurements alone is considered in our decorrelation approach.

First, we constructed L matrix to decorrelate the residual vector, with respect to its covariance matrix. The residual gain matrix is now equivalent to the weight matrix used to decorrelate the residual vector.

The covariance matrix of a residual vector can be computed according to Equation (12) as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16)

The residual vector is decorrelated using its covariance matrix, because correlation among residual vectors will be significant, when datasets contain many outliers. The covariance matrix of a residual vector can be determined from the observation weight matrix and the covariance matrix of the innovation vectors; therefore, the residual vector of covariance is relevant. Conversely, we can determine the correlation of the residual vectors based on the residual covariance correlation. First, we constructed a matrix L to meet the condition det (L) = 1, and then decorrelated the covariance matrix of the residual vector, as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (17)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is a diagonal matrix. This equation can also be solved by adopting the Cholesky decomposition.

Similarly, we applied the L decorrelated matrix of the residual vector to obtain the decorrelated residual vector [V.sub.decorrelation], as follows:

[V.sub.decorrelation]=LV (18)

Finally, an iterative algorithm was used to reduce the residual correlation of the residual vector, according to the following steps:

1. Using Equation (17) and the covariance matrix of the residual vector, the decorrelated matrix L is solved;

2. Using the L matrix and the residual vector, the decorrelated residual vector (Equation (18)) is calculated;

3. Using the decorrelated residual vector (Equation (15)), robust classification factors is used to calculate a robust factor. Subsequently, the equivalent weight matrix is computed (Equations (13) and (14));

4. The gain matrix is calculated by the equivalent weight matrix (Equation (9));

5. The residual vector based on this gain matrix is solved (Equation (11));

6. Steps (3) to (5) are repeated, until the residual vector correction value is less than a given threshold value.

This process is illustrated as a schematic in Figure 1.

4. EXPERIMENTAL ANALYSIS

The observational data and precision products provided by IGS station BJFS in China were used to evaluate the proposed method in this study. Specifically, the observations of 1,000 epochs on August 29, 2010 during the period of 00:00-08:20 comprised our sample dataset. The sampling interval was 30 s. The MW/TECR method (Liu, 2011) was used to detect and repair the cycle slip in these original observations. Pre-processing of these data indicated that few inherent outliers occurred before 300 epochs; hence, these data were used as simulation data to verify our proposed algorithm. The data in the range of 300 to 1000 epochs contained many divergent outliers, which were identified using the traditional Kalman filter algorithm.

To test the validity of our algorithm, the coordinates of this IGS station at the SOPAC site was taken as the true values (http://sopac.ucsd.edu/sector.shtml). Moreover, the true value was compared with the coordinates determined using different algorithms. The Minimally Detectable Bias(MDB) (Kouba and Heroux, 2001) used by many scholars to determine outlier size gave values of 3 cm and 4 m for the carrier phase and pseudorange observations in this dataset, respectively.

As the first step, we added outliers of a certain size to the PRN14 satellite data during the 100th epoch to determine the influence of outliers on the final positioning result. Figure 2 shows both the innovation and the residual values for the raw data of the carrier phase and the pseudorange with added outliers of 0.1 and 10 m for the signal satellite in the 100th epoch, respectively. The innovation and the residuals vectors are subject to Gaussian distributions. In addition, the innovation vector of the carrier phase is defined at the decimeter to meter-level, whereas the residual is defined at the centimeter-level. In contrast, the innovation vectors of pseudorange and its residuals are defined at the meter-level. The innovation values of all satellites in the carrier phase are similar within the same epoch; however, the innovation vectors and residuals of satellites with outliers are abnormal, when the carrier phase contains 0.1 m-outliers. We can determine the location and size of these outliers based on the innovation value of the carrier phase containing these outliers. However, this innovation vector cannot be derived accurately when many outliers are detected. Unlike the innovation vector, residuals can be clearly detected, even when the observation data of a satellite contains outliers, as shown in Figure 2a (2). Finally, the residual vectors of all satellites show marked changes (Fig. 2a (2)), when 0.1 m small outliers occur in the carrier phase observation values of a single satellite. This result suggests that all satellite residuals are correlated to a certain extent.

Figure 3 depicts the robust factor of the innovation and the residual vectors. The innovation values of all satellite carriers are similar within the same epoch, when the carrier phase observation value is very precise. In this case, we cannot separate outliers from correct observations effectively when there are many outliers present; hence, the robust factor of the innovation vector cannot resist the effects of outliers. Given the low-accuracy of pseudorange observations, both innovation vectors and residual vectors can still be detected effectively, and the number of satellites with outliers can be readily identified. The process of determining correlations among the robust factors of the residual vector is relatively simple (Fig. 3a (2)); therefore, we are able to eliminate the influence of this vector and its robust factor.

Figsures 4a and 4b display the residual vector and its robust factor after decorrelation using Equations (17) and (18). The correlation between the residual vector and its robust factor is greatly reduced using the L transformation, showing that the decorrelation between the L matrix and the residual vector is more effective than other correlations. In contrast, the L matrix based on the covariance of the residual vector is not very accurate. Thus, the method shown in Figure 1 will greatly weaken the residual vector and its robust factor.

Figsures 5a and 5b exhibit the residual vector and its robust factor, after the decorrelation and iteration proposed in this study. The correlation between the residuals of each satellite in the carrier phase and its robust factor is eliminated.

To demonstrate the effectiveness of our algorithm, different schemes were used to explore outliers of a specific size in the 100th epoch; we carried out the following modifications to the original dataset:

1. The carrier phase had 0.1 m-outliers added to its values;

2. The pseudorange had 10 m-outliers added to its values;

3. The carrier phase observations had additional 0.5 m-outliers, and the pseudorange had 50 moutliers added to their values.

The following algorithms were used to compare and analyze the original and modified datasets that with or without outliers:

Scheme 1: Data were processed using a standard Kalman filter.

Scheme 2: Data were processed using the RKF based on innovation vectors.

Scheme 3: Data were processed using the RKF based on residuals, involving decorrelation and iteration for the residual vector and its robust factor.

Table 1 presents the results of these three schemes, showing the deviation of the estimated coordinate, given outliers to the 100th epoch of the raw data. The standard Kalman filter was clearly ineffective against outliers, and the effect of the 0.1 m-outliers in the carrier phase observations on the position parameters was much stronger than the 10 m-outliers in the pseudorange observations. Scheme 2 dealt robustly for the pseudorange observations with outliers, but showed large deviations for the carrier phase observations with outliers. The innovation vector was insensitive to the carrier phase outliers, suggesting that scheme 2 is ineffective against small outliers occurring in highly accurate observations. The coordinate deviations in the scheme 3 were generally small (Fig. 5), suggesting that the correlation among the carrier phase observation values of various satellites is greatly reduced by this scheme. In theory, the accuracy of scheme 3 is higher than scheme 2, although scheme 2 is better at addressing outliers.

Our results show that the coordinate deviations resulting from introduced outliers are mainly influenced by carrier phase observations. Figure 6 shows the residual time series of 3D coordinates obtained from the three schemes described above in the 100th epoch, when 0.1 m-outliers were introduced in the original carrier phase observation data. The 3D coordinate residuals of schemes 1 and 2 show a marked change in the 100th epoch. In contrast, scheme 3 did not undergo such a large change at this point, showing that it can effectively resist the effect of such outliers. Thus, scheme 3 can effectively resist the influence of a single outlier in a single satellite observation.

Above content, we analyzed and compared the different schemes applied to a single satellite observation set. The actual measured data contained many outliers within the same epoch; different satellites may also contain outliers of various sizes. The following three schemes were used to investigate the PPP problem, with the solution to the coordinate time series shown in Figure 7. Figures 7a, 7b, and 7c show X, Y, and Z directions, respectively.

Scheme 4: Data processed using the standard Kalman filter.

Scheme 5: Data processed using the RKF based on residual, with decorrelation for the residual and its robust factor.

Scheme 6: Data processed using the RKF based on residuals, with decorrelation, and iterations for its residual and robust factor.

Scheme 5 can effectively resist the few outliers in observation values occurring before 300 epochs, because this scheme reduces the weight of most of the observed values. In contrast, the coordinate solution has little effect after 300 epochs. The RKF does not appear to diverge, either when the number of outliers increases, or when multiple satellites have outliers within the same epoch. Nonetheless, scheme 6 is more effective than scheme 5 against these outliers. In spite of the abundance of outliers, the filters did not diverge, mainly because the correlation between residual vectors of different satellites in the carrier phase was effectively reduced. In view of this limitation, distinguishing outliers from correct observation values based on the innovation vector alone is difficult, while the residual vector can be used to determine whether or not a large difference exists because of outliers. In addition, the coordinate series for scheme 6 clearly does not diverge (Fig. 7). Furthermore, its convergence rate is the fastest, and coordinates of this scheme deviate least from the true coordinates.

Table 2 compares the accuracies of schemes 4, 5, and 6 from the 300th to 1000th epoch. The scheme 4 does not show any effect of the difference among outliers; moreover, the scheme 5 can be more effective to resist outliers. Our experimental analyses show that using our proposed method (scheme 6) effectively resists various sizes of outliers in various positions in both carrier phase and pseudorange observations, and have the minimal STD and RMS errors.

5. CONCLUSIONS

An RKF method based on the residual vector is proposed for static PPP in this study. Experimental analyses show that outliers can be identified using the innovation vector, containing only one observation. Nonetheless, outliers are difficult to distinguish in unequal precision observations using innovation vectors. In contrast, the residual vector effectively reflects outliers in such observations, because of its dependence on the state parameters of the posteriori case. Furthermore, residuals are more sensitive to outliers than innovation vectors; the robust factor of these residuals resists the influence of outliers on PPP effectively.

The residual component is affected by any single observation outliers; hence, decorrelation and iteration must be carried out. This decorrelated residual vector accurately reflects the position and size of outliers in observations. Our experimental results show that the proposed RKF method resists outliers more effectively and accurately than any other traditional RKF method with unequal precision observations.

ACKNOWLEDGEMENT

This work has been supported by the Fundamental Research Funds for the Central Universities (grant number: 2014ZDPY29) and a Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions (grant number: SZBF2011-6-B35). The author is grateful for the IGS community for providing GNSS data and precision products in this study. The anonymous reviewers and the Editor Dr. Zdenka Schenkova are thanked sincerely for their constructive comments and suggestions to improve the quality of the paper.

REFERENCES

Chang, G.: 2014a, Kalman filter with both adaptivity and robustness. Journal of Process Control, 24, No. 3, 81-87. DOI: 10.1016/j.jprocont.2013.12.017

Chang, G.: 2014b, Robust Kalman filtering based on Mahalanobis distance as outlier judging criterion. Journal of Geodesy, 88, No.4, 391-401. DOI: 10.1007/s00190-013-0690-8

Chang, G. and Liu, M.: 2015, An adaptive fading Kalman filter based on Mahalanobis distance. Proceedings of the Institution of Mechanical Engineers Part G: Journal of Aerospace Engineering, 229, No.6, 1114-1123. DOI: 10.1177/0954410014545181

Collins, P., Bisnath, S., Lahaye, F. and Heroux, P.: 2010, Undifferenced GPS ambiguity resolution using the decoupled clock model and ambiguity datum fixing. Navigation, 57, No .2, 123-135. DOI: 10.1002/j.2161-4296.2010.tb01772.x

Dawidowicz, K. and Krzan, G.: 2014, Coordinate estimation accuracy of static precise point positioning using online PPP service, a case study. Acta Geodaetica et Geophysica, 49, 37-55. DOI: 10.1007/s40328-013-0038-0

Hefty, J. and Gerhatova, L.: 2012, Potential of precise point positioning using 1 Hz GPS data for detection of seismic-related displacements. Acta Geodyn. Geomater., 9, No. 3 (167), 303-313.

Hu, H., Gao, J. and Yao, Y.: 2014, Land deformation monitoring in mining area with PPP-AR. International Journal of Mining Science and Technology, 24, No. 2, 207-212. DOI: 10.1016/j.ijmst.2014.01.011

Huber, P.J.: 2011, Robust statistics. Springer Berlin Heidelberg.

Gao, Y. and Shen, X.: 2002, A new method for carrier phase based on Precise Point Positioning. Navigation, 49, No. 2, 109-116. DOI: 10.1002/j.2161-4296.2002.tb00260.x

Ge, M., Gendt, G., Rothacher, M., Shi, C. and Liu, J.: 2008, Resolution of GPS carrier-phase ambiguities in Precise Point Positioning (PPP) with daily observations. Journal of Geodesy, 82, No.7, 389-399. DOI: 10.1007/s00190-007-0187-4

Guo, F. and Zhang, X.: 2014, Adaptive robust Kalman filtering for precise point positioning measurement. Science and Technology, 25, 105011-105018. DOI: 10.1088/0957-0233/25/10/105011

Kalman, R.E.: 1960, A new approach to linear filtering and prediction problems. Journal of Fluids Engineering, 82, No. 1, 35-45. DOI: 10.1115/1.3662552

Kouba, J. and Heroux, P.: 2001, Precise point positioning using IGS orbit and clock products. GPS solutions, 5, No. 2, 12-28. DOI: 10.1007/PL00012883

Krzan, G. and Przestrzelski, P.: 2016, GPS/GLONASS precise point positioning with IGS real-time service products. Acta Geodynamica et Geomaterialia, 13, No. 1(181), 69-81. DOI: 10.13168/AGG.2015.0047

Laurichesse, D., Mercier, F., Berthias, J.P., Broca, L. and Cerri, L.: 2009, Integer ambiguity resolution on undifferenced GPS phase measurements and its application to PPP and satellite precise orbit determination. Navigation, 56, No. 2, 12-28. DOI: 10.1002/j.2161-4296.2009.tb01750.x

Li, X., Ge, M., Dai, X., Ren, X., Fritsche, M., Wickert, J. and Schuh, H.: 2015a, Accuracy and reliability of multi-GNSS real-time precise positioning: GPS, GLONASS, BeiDou, and Galileo. Journal of Geodesy, 89, No. 6, 607-635. DOI: 10.1007/s00190-015-0802-8

Li, X., Ge, M., Dousa, J. and Wickert, J.: 2014, Real-time precise point positioning regional augmentation for large GPS reference networks. GPS Solutions, 18, No. 1, 61-71. DOI: 10.1007/s10291-013-0310-3

Li, X., Zhang, X., Ren, X., Fritsche, M., Wickert, J. and Schuh, H.: 2015b, Precise positioning with current multi-constellation Global Navigation Satellite Systems: GPS, GLONASS, Galileo and BeiDou. Scientific reports, 5, 1-14. DOI: 10.1038/srep08328

Liu, Z.: 2011, A new automated cycle slip detection and repair method for a single dual-frequency GPS receiver. Journal of Geodesy, 85, No. 3, 171-183. DOI: 10.1007/s00190-010-0426-y

Nie, J., Zhang, S., Xu, Y., Zhang, Y. and Wang, Y.: 2010, Precise point positioning based on robust Kalman filtering. Journal of Earth Sciences and Environment, 32, No. 2: 218-220.

Shi, J. and Gao, Y.: 2014, A comparison of three PPP integer ambiguity resolution methods. GPS Solutions, 18, No. 4, 519-528. DOI: 10.1007/s10291-013-0348-2

Trasak, P. and Stroner, M.: 2014, Outlier detection efficiency in the high precision geodetic network adjustment. Acta Geodaetica et Geophysica, 49, No. 2, 161-175. DOI: 10.1007/s40328-014-0045-9

Xu, C., Gao, J., Hu, H. and Wang, J.: 2012, Robust Kalman filtering for Precise Point Positioning. Journal of China University of Mining & Technology, 41, No. 5, 857-862. DOI: 10.1109/JSEE.2013.00081

Yang, Y.: 1997, Robust Kalman filter for dynamic systems. Joumal of the PLA Institute of Surveing and Mapping. 13, No. 2, 79-84.

Yang, Y. and He, H.: 2001, Adaptive robust filtering for kinematic GPS positioning. Acta Geodaetica et Cartographica Sinica, 30, No. 4, 293-298.

Yang, Y., Song, L. and Xu, T.: 2002, Robust estimator for correlated observations based on bifactor equivalent weights. Journal of Geodesy, 76, No. 6, 353-358. DOI: 10.1007/s00190-002-0256-7

Yang, Y., He, H. and Xu, G.: 2001, Adaptively robust filtering for kinematic geodetic positioning. Journal of Geodesy, 75, No. 2-3, 109-116. DOI: 10.1007/s001900000157

Zhang, X., Guo, F. and Zhou, P.: 2014, Improved precise point positioning in the presence of ionospheric scintillation. GPS Solutions, 18, No. 1, 51-60. DOI: 10.1007/s10291-012-0309-1

Zhang, X., Guo, F., Li, P. and Zuo, X.: 2012, Real-time quality control procedure for GNSS precise point positioning. Geomatics and Information Science of Wuhan University, 37, No. 8, 940-944.

Zumberge, J.F., Heflin, M.B., Jefferson, D.C., Watkins, M.M. and Webb, F.H.: 1997, Precise point positioning for the efficient and robust analysis of GPS data from large networks. Journal of Geophysical Research: Solid Earth, 102, No. B3, 5005-5017. DOI: 10.1029/96JB03860

Wang, J., Wang, J. and Gao, J.: 2008, Study on GNSS navigation model based on EKF. Zhongguo Kuangye Daxue Xuebao/Journal of China University of Mining and Technology, 37, No. 4, 473-477.

Yi-fei YAO (1), Jing-xiang GAO (1)(*), Zeng-ke LI (1), Chang-hui XU (2) and Xin-yun CAO (3)

(1) NASG Key Laboratory of Land Environment and Disaster Monitoring, China University of Mining & Technology, Xuzhou, Jiangsu 221116, China

(2) Chinese Academy of Surveying & Mapping, Beijing 100830, China

(3) Wuhan University, Wuhan, Hubei 430072, China

Cite this article as: Yao Y, Gao J, Li Z, Xu Ch, Cao X: A new robust Kalman filtering algorithm of unequal precision observations based on residual vectors in static precise point positioning. Acta Geodyn. Geomater., 13, No. 4 (184), 397-408, 2016. DOI: 10.13168/AGG.2016.0022

(*) Corresponding author's e-mail: jxgao@cumt.edu.cn

ARTICLE INFO

Article history:

Received 31 March 2016

Accepted 23 June 2016

Available online 13 July 2016

Keywords:

Robust Kalman Filter (RKF)

Precise Point Positioning (PPP)

Innovation vector

Residual vector

Unequal precision

Carrier phase observations

Table 1 Coordinate deviations resulting from different schemes on observations having outliers. Schemes The size of added outliers The coordinate deviation X/mm Y/mm Z/mm Carrier phase 0.1m; 0 37 46 Schemes 1 Pseudorange 10m; 3 2 2 Carrier phase 1m, pseudorange 29 395 484 100m; Carrier phase 0.1m; 0 19 24 Schemes 2 Pseudorange 10m; 0 0 0 Carrier phase 1m, pseudorange 100m; Carrier phase 0.1m; 0 0 -1 Schemes 3 Pseudorange 10m; 0 -1 0 Carrier phase 1m, pseudorange 100m; Table 2 Comparison of the precision and accuracy of schemes 4, 5 and 6. scheme STD RMS X/mm Y/mm Z/mm X/mm Y/mm Z/mm scheme 4 81 112 93 121 111 82 scheme 5 7 20 16 6 91 93 scheme 6 10 15 16 12 26 62

Printer friendly Cite/link Email Feedback | |

Title Annotation: | ORIGINAL PAPER |
---|---|

Author: | Yao, Yi-Fei; Gao, Jing-Xiang; Li, Zeng-Ke; Xu, Chang-Hui; Cao, Xin-Yun |

Publication: | Acta Geodynamica et Geromaterialia |

Article Type: | Report |

Date: | Oct 1, 2016 |

Words: | 6632 |

Previous Article: | The influence of different types of noise on the velocity uncertainties in GPS time series analysis. |

Next Article: | Radon release from underground strata to the surface and uniaxial compressive test of rock samples. |

Topics: |