2. AMBIGUITY RESOLUTION AND VALIDATION

2.1. Ambiguity Resolution

Carrier phase ambiguity resolution is the key to high accuracy positioning with GNSS. Reliable ambiguity resolution is a function of several factors, the main ones being type of measurement, residual measurement errors, geometry and algorithm formulation. To have a good chance of reliable ambiguity resolution, the overall error budget in the observation is generally required to be at the half a cycle level with an uncertainty of less than a quarter of a cycle (Sauer, Reference Sauer2003). This requirement can be achieved through the use of:

• • Linear combinations of raw measurements to form longer wavelengths to aid the resolution of ambiguities for shorter wavelengths (e.g. the Three-Carrier Ambiguity Resolution – TCAR).

• • Single/double differencing to remove common errors and mitigate correlated errors.

• • External products to mitigate errors (e.g. satellite orbits and clocks provided by the International GNSS Service [IGS]).

• • Other error corrections from dedicated networks (e.g. Un-calibrated Phase Delays – UPD).

• • A combination of the above.

After pre-processing, the model used for GNSS positioning is:

(1)

where y is the observation vector, a and b are unknown parameter vectors, e is the noise vector, a∈Zⁿ is the integer ambiguity, contains the other parameters to be solved (e.g. position), and .

The estimation of integer values of a is not straightforward. The first step in the Integer Least Squares (ILS) method (Teunissen, Reference Teunissen1993) is to use the traditional least squares approach to estimate a as real (float) values, â. Denoting the corresponding solution of b as , the vector of observation residuals can be written as:

(2)

The Sum of the Squared Error (SSE) is given by:

(3)

The matrix G _y is the cofactor matrix of the variance-covariance matrix Q _y.

If ê∈N(0, σ ²), then

(4)

(5)

The next step is to map the ambiguity from an n-dimensional real space to n-dimensional integer space. There are various methods for the mapping, with the simplest being rounding to the nearest integer. When a number of carrier phase observations are involved, a float vector of ambiguities can be used as the initial vector. The corresponding integer vector can then be obtained by rounding each element of the float vector to its nearest integer. This method has the disadvantage that it does not take into account any of the correlation that may exist between the individual elements of the ambiguity vector. In this sense, this simple method can be only safely used if the model is improved by using additional information such as precise ionospheric corrections (Hernández-Pajares et al., Reference Hernández-Pajares, Juan, Sanz, Aragón, Ramos, Samson, Tossaint, Albertazzi, Odijk, Teunissen, de Bakker, Verhagen and van der Marel2010).

Another relatively easy way to obtain an integer ambiguity vector from the float ambiguity vector is to apply a sequential rounding scheme to the elements of the latter. This approach uses rounding to determine the integer ambiguity for the first element of the float vector, which typically has the smallest estimated error, and therefore, the best chance for the correct integer ambiguity. The remaining ambiguities are then sequentially rounded to the nearest integers after taking into account their correlations with the others that have been resolved. This is referred to as integer bootstrapping. The advantage of this method is its simplicity. However, its results depend on the parameterisation of the ambiguities. This implies that results will differ when applying this estimator to either original or decorrelated ambiguities.

Another approach for estimating an integer ambiguity vector using the integer least squares estimator is by solving the following minimization problem (Teunissen, Reference Teunissen1993):

(6)

This estimator is optimal in the sense that it gives the highest probability of finding the correct integer vector. The estimator can be interpreted as finding the integer vector that has the shortest distance to the float solution, measured in the metric of the variance-covariance matrix Q _â. This method does not have the disadvantages of rounding and bootstrapping. However, it is more complex. Moreover, a solution using standard least-squares algorithms cannot be found because of the integer nature of the solution. A discrete search through the complete space of integers Z _n may be required in order to obtain optimal results. The different elements of the float ambiguity vector may be highly correlated. Hence, for methods such as rounding and bootstrapping, the more correlated the ambiguities are, the more likely that these methods will yield non-optimal results.

A brute force search is a time consuming procedure and may not be acceptable for some applications including those where time is critical or where cost is a concern. A procedure to reduce the search space developed by Teunissen (Reference Teunissen1993), namely decorrelation of the ambiguities, is widely accepted by the GNSS community. Based on the least squares ambiguity estimation, a transformation of the ambiguities and corresponding variance-covariance matrix into an equivalent but less correlated set of ambiguities is carried out before estimating the ambiguities as integers. The so-called Z transformation is expressed as:

(7)

(8)

Z is a matrix with all its elements as integers, and the inverse of the Z-matrix should exist. The Z transformation should provide maximum decorrelation of the ambiguities.

The more recent Integer Aperture (IA) estimation method defines the size and shape of the aperture pull-in regions (Teunissen, Reference Teunissen2003, Reference Teunissen2005). In classical hypothesis testing theory, the size of an acceptance region is determined by the choice of the testing parameters: the false alarm rate and missed detection probability (or detection power). However, in the case of integer ambiguity resolution, the selection of these parameters is not obvious. It is especially important that the probability of incorrect ambiguity fixing is small. Therefore, the concept of IA estimation with a fixed fail rate has been introduced. This means that the size of the aperture space is determined by the condition that the fail rate is either equal to or lower than a fixed value. At the same time, the shape of the aperture pull-in regions should preferably be chosen such that the success rate is still as high as possible.

After the process of optimisation employing the techniques above, the estimate of a as integer values is denoted as and the corresponding solution of b is denoted as . The vector of observation residuals can then be written as _, where i is the order of a candidate ambiguity. The number of the candidates depends on the search space. However, only the best and second best set of ambiguities are considered for validation. Therefore, 0<i≤2. This is explained further in the next section.

The Sum of the Squared Error (SSE), , can then be expressed as:

(9)

If , then

(10)

(11)

After the determination of the ambiguity candidates, the ambiguity residuals are determined as:

(12)

The SSE, R _i, is expressed as:

(13)

and,

(14)

where δ _i is the non-central parameter of Chi-square distribution.

2.2. Validation of the Ambiguities

Each integer estimates of the ambiguities contain a degree of uncertainty for both integer rounding (e.g. TCAR) and least squares methods (e.g. LAMBDA). Therefore, the resolved ambiguities must be validated before they can be used for high accuracy positioning. In general, validation is based on the formation of a test statistic and the determination (or assumption) of its distribution. This distribution should in theory enable the threshold or confidence level to be computed for comparison with the test statistic, in order to determine if the integer ambiguity is acceptable. There are a number of tests based on the SSE of either the observation or ambiguity residuals (the latter is a subset of the former).

2.2.1. Validation Based on Observation Residuals

The ratio test (also referred to as the F-ratio test) constructed based on the observation (post-fit) residuals between the second best and best is

(15)

where, k is the threshold. The best and second best can also be tested separately as (Chen, Reference Chen1997):

(16)

where χ _γ²(m−p, 0) is the critical value corresponding to the central Chi-square distribution with a level of significance γ (probability in the right-hand tail); χ _β²(m−p, δ) is the critical value corresponding to the non-central Chi-square distribution with a probability β and non-centrality parameter:

(17)

The difficulty in using this test is the choice of the critical values of γ and β.

As an alternative, the W-ratio test was proposed by Wang et al. (Reference Wang, Stewart and Tsakiri2000):

(18)

where and var(d)=δ ²Q _d.

Q _d is the variance co-factor of d and δ ² is the so-called variance factor. However, this approach makes the incorrect assumption that the use of a posteriori variance cofactor translates into the W-ratio with a Student's t distribution (Verhagen, Reference Verhagen2004).

2.2.2. Validation Based on Ambiguity Residuals

Similar to the F-ratio test, Euler and Schaffrin (Reference Euler and Schaffrin1991) proposed a test given by:

(19)

where, k is the threshold. Notably, a number of constant values for k have been used without a credible theoretical or practical justification. For example, the values of 1·5 and 3 have been used by Wei and Schwarz (Reference Wei and Schwarz1995), and Leick (Reference Leick2003) respectively.

Another approach is based on the difference of the quadratic forms R ₁ and R ₂ (Tiberius and De Jonge, Reference Tiberius and De1995). Firstly, the test is carried out. When this is passed, the next step is to perform the difference test:

(20)

The use of tests based on either ratio or difference involving the best and second best ambiguity vector, enables all cases to be covered, as they represent the smallest values (i.e. deviations). However, the difficulty in using this test lies in the choice of the critical value, k which is determined empirically.

A class of Integer Aperture (IA) estimation and validation methods is proposed by Teunissen (Reference Teunissen2003, Reference Teunissen2005). It is very important to note that IA estimation with a fixed fail rate is a method that involves both integer estimation and validation, and allows for an exact and overall probabilistic evaluation of the solution. With the traditional approaches (e.g. the conventional ratio test applied with a fixed critical value or threshold) an overall probabilistic evaluation of solution is not possible.

If the float ambiguity lies in one of the acceptance regions, the corresponding resolved integer solution is accepted. The size of the regions depends on the choice of the threshold value μ (i.e. the larger it is, the larger the acceptance region becomes). In the limiting case, μ=1, the acceptance regions are equal to the integer least-squares pull-in regions, and hence, the integer solution is always accepted. In the other limiting case, μ=0, the integer solution is always rejected.

IA validation is based on a fixed success/fail rate. The method uses the reciprocal of expression (19) as the test statistic. It involves the derivation of the reciprocal of expression (19), μ ₀ with measurements at the current epoch, and using Monte Carlo simulation to generate float ambiguities and corresponding reciprocals of expression (19), μ _i (Teunissen and Verhagen, Reference Teunissen and Verhagen2004). Therefore, a statistic for the success/fail rate can be calculated using μ ₀ and μ _i together with other parameters. In order to obtain a good approximation, a Monte Carlo simulation with a large number of samples is required (typically >100,000). This is a major limitation for this method for real time applications.

From the preceding sections, it is clear that ambiguity validation is an open problem. This is because none of the available integer validation test statistics is based on sound theoretical principles. Furthermore, there is no single test that can be used in all situations (Verhagen, Reference Verhagen2004). A common problem with the observation and ambiguity residuals based tests above, is the choice of the threshold. Either an empirically determined value is used, or an assumption is made on the distribution. Although the IA estimator based validation is a promising approach, it has some weaknesses in practice.

3. NEW APPROACH TO CARRIER PHASE BASED INTEGRITY MONITORING

In order to detect failures early, the approach to carrier phase based integrity monitoring proposed in this paper involves two steps executed at the ambiguity resolution and positioning stages (the latter comprising the detection function in the measurement domain and the use of the protection levels in the position domain). This section addresses the details of the integrity monitoring at the ambiguity resolution stage. The second stage is based on the Carrier phase RAIM (CRAIM) developed for the conventional RTK (Feng et al., Reference Feng, Ochieng, Moore, Hill and Hide2009). This is easily transferable to other positioning concepts including PPP. This is elaborated in the next section.

The ambiguity validation process can be considered as a failure detection process in the ambiguity domain. Four stages are normally involved in the process: the construction of test statistics, description of the distribution of the test statistics, determination of threshold, and determination of integrity flag. The test statistics formulated are functions of a number of factors including geometry, observations, observation residual errors after pre-processing, and noise. The threshold should reflect the correctness/confidence level (or success rate) of the test in order that decisions are made based on it.

As discussed in the preceding section, ratio based testing is commonly used for ambiguity validation. The test statistic usually takes the form of either expression (19) or its reciprocal and a fixed threshold. Furthermore, the fixed success/fail rate based validation approach employed with the IA method has limitations associated with real-time applications. Therefore, the remaining major issue is the determination of threshold or confidence level which is based on an accurate description of the distribution of the test statistic to be used instead of the fixed thresholds used currently.

The approach in this paper uses the same test statistic as that used in conventional ratio test (expression 19), but derives a distribution for the test statistic (discussed below) and uses it within a numerical algorithm to compute the confidence level.

In order to determine the confidence level of the ratio test, the distribution of (R ₂/R ₁) is required. Thus expression (19) can be rewritten as:

(21)

As can be seen from expression (14), both the numerator and denominator follow a non-central χ ² distribution. Therefore, if R ₁ and R ₂ are independent (assumed in this paper) then (R ₂/R ₁) has a doubly non-central F-distribution (Bulgren, Reference Bulgren1971).

(22)

Therefore, the confidence level p _c of (R ₂/R ₁) can be derived from:

(23)

where n is the number of ambiguities; δ ₁ and δ ₂ can be determined using the traditional failure detection scheme with a probability of false alert (P_FA) and a probability of missed detection (P_MD).

In expression (14), there is one P_FA and two probabilities of missed detection (P_MD1 and P_MD2) for the best and second best SSE, R _i(i=1, 2). δ ₁ and δ ₂ can then be derived from the following expression.

(24)

The analytical method for the solution of δ ₁ and δ ₂ is not straightforward, requiring the application of numerical methods using series representation. This is also used to determine the confidence level in expression (23). The CDF of the general expression (23) can be re-written as (Bulgren, Reference Bulgren1971)

(25)

where

• • and are the probabilities of Poisson distribution, Γ( ) is the Gamma function,

• • is the CDF of Beta distribution with u=n ₂x/(n ₂x+n ₁) and x⩾0, is the Beta function with a>0 and b>0.

For a feasible implementation of the algorithm, the two infinite series are truncated when the higher order is not needed in reaching a specified accuracy.

Following the execution of LAMBA type algorithms, the method here takes the outputs of carrier phase ambiguity resolution (R ₁, R ₂ and n) together with pre-defined P_FA, and P_MDs, to generate in an online computation, the confidence levels for the candidate ambiguities at each epoch. Given a threshold, the confidence levels can be used to accept or reject the best set of candidate ambiguities. Some users may be more concerned with the confidence level satisfying the requirements rather than the confidence level itself. In this case, in order to reduce the processing resources consumed by online computation, it is quite possible to carry out offline calculation for the threshold of (R ₂/R ₁) to reflect the required confidence. Therefore, a look-up table can be calculated offline taking all possible cases including selecting a range of confidence levels and degrees-of-freedom. This is especially useful for embedded systems where computational power is low.

4. APPLICATION OF NEW APPROACH TO PPP

Crucial to PPP is the requirement to deal with all the error sources including those that are either eliminated or reduced by differenced observables in conventional RTK. Therefore, the PPP user algorithm is designed, in this paper, to include specific error models including Earth tides and wind up. The algorithm is based on sequential least squares that exploits the current and previous data. Figure 1 shows a high level architecture of the PPP user algorithms that employs the two-step integrity monitoring process (ambiguity resolution and positioning stages) highlighted in yellow and green respectively. The first stage employs the new validation algorithms presented in the previous section. The second stage is elaborated in this section as it applies to PPP.

Figure 1. High level architecture of PPP and integrity monitoring.

In a filtering based carrier phase positioning, the ambiguity vector is considered as unknown in the state equation. If the ambiguity vector is resolved correctly, the integer ambiguity vector therefore substitutes the float values in the state. The ambiguity resolution and validation can be carried out separately with the input from the filter. Therefore, the integrity monitoring method proposed and presented in the previous section is executed at this stage.

For the integrity monitoring at the positioning stage, the conventional RTK based method (CRAIM) can be migrated to PPP. The information needed include those items that can be extracted from the PPP Kalman filter, given by product providers and agreed values by users on requirements and the other residual errors. The outputs from the integrity monitoring algorithm are the protection levels and integrity status.

The information that can be extracted from the PPP Kalman filter are (Feng et al., Reference Feng, Ochieng, Moore, Hill and Hide2009):

• • The innovation vector r

• • The weight matrix W

• • The design matrix H

• • The measurement noise matrix R

• • The gain K

• • The covariance matrix P

The information given by product providers are clock and orbit corrections for each satellite and associated residual uncertainty σ _i,SIS, and corrections for the ionospheric error and the associated residual uncertainty σ _i,IONO. If the residual uncertainty information is not available, a value agreed by users which overbounds these residual errors can be used. In addition, σ _i,User is the standard deviation of residual errors for each measurement after applying appropriate models and includes residual multipath, and receiver noise. The values of σ _i,User which overbound these residuals have to be agreed by users. Furthermore, the values of the parameters of the Required Navigation Performance (RNP) have to be provided by users, especially the integrity and continuity parameters used in this part.

In PPP (Figure 1), the observables used include ionospheric free pesudorange (P_C), ionospheric free carrier phase (L_C), Melbourne-Wubbena (L_W-P_n), and geometry free linear combinations (LI). These observables reduce the impact of bias type error on PPP performance, although the noise levels are amplified. The standard deviation of each observable is derived from error propagation theory (σ _obs=f (σ _obs,user, σ_SIS, σ_{obs, IONO})).

The methods for determining the test statistics, thresholds and protection levels are similar to those used in CRAIM (Feng et al., Reference Feng, Ochieng, Moore, Hill and Hide2009). The construction of the test statistics consisting of one full set and a number of subsets is based on the parameters extracted from the PPP Kalman filter. The test statistics derived from the innovation of a Kalman filter follow a normalized Chi distribution with degrees of freedom equal to the number of measurements used. Whether it is a central or non-central Chi-squared distribution depends on the absence or presence of a range bias error (Diesel, Reference Diesel and Luu1995; Lee, Reference Lee and Laughlin1999; Feng et al., Reference Feng, Ochieng, Walsh and Ioannides2006). Therefore, the corresponding thresholds can be determined from the Chi distribution for a given probability of false alert. To detect the presence of failure, the test statistics are compared to corresponding thresholds. If there is a failure, either an integrity flag is raised or a failure exclusion process is invoked. It should be noted that the innovation sequence contains information obtained from the previous states. Therefore, the response to failure will be delayed for a few epochs and the same for positioning accuracy.

There are two sets of protection levels determined in the PPP integrity monitoring algorithm. One set is determined by using the horizontal and vertical uncertainty in the positioning solution which can be extracted from the covariance matrix of the PPP Kalman filter. However, the horizontal and vertical uncertainties do not immediately reflect the cases when swapping between float and fixed ambiguities. This is because of sudden changes in the positioning accuracy. A few epochs may be required for the covariance matrix to adjust. This needs to be considered in the overall protection level determination. The other set of protection levels is based on the slope concept. Similar to conventional RAIM, the second set of protection levels has little relation to the residuals. In the case that two sets of protection levels are available, the larger ones are adopted. Both the horizontal and vertical protection levels are used to compare against alert limits specified by the user. In safety critical applications, if the protection level is larger than the alert limit, a warning (alert) should be issued to let the user know of the integrity status. However, in availability critical applications, the alert limit may be relaxed so that the user can continue to use the system without being interrupted by an integrity alert.