Single-receiver single-channel multi-frequency GNSS integrity: outliers, slips, and ionospheric disturbances

In this contribution the integrity of single- receiver, single-channel, multi-frequency GNSS models is studied. The uniformly most powerful invariant test statistics for spikes and slips are derived and their detection capabilities are described by means of minimal detectable biases (MDBs). Analytical closed-form expressions of the phase-slip, code-outlier and ionospheric-disturbance MDBs are given, thus providing insight into the various factors that contribute to the detection capabilities of the various test statistics. This is also done for the phaseless and codeless cases, as well as for the case of a temporary loss-of-lock on all frequencies. The analytical analysis presented is supported by means of numerical results.

multi-receiver/antenna case, sometimes even with additional constraints included. An example of the latter is the quality control of baseline-constrained GNSS attitude models (Giorgi et al. 2012), while geometry-dependent receiver autonomous integrity monitoring (RAIM) is an example of the single-receiver, multi-channel case (Teunissen 1997;Wieser et al. 2004;Feng et al. 2006;Hewitson and Wang 2006).
In the present contribution, we study the single-receiver, single-channel model. It is the weakest model of all, due to the absence of the relative receiver-satellite geometry. Despite its potential weakness, there are several advantages to single-receiver, single-channel data validation. First, since it is the simplest model of all, it can be executed in real-time inside the (stationary or moving) receiver, thus enabling early quality control on the raw data. Second, the geometry-free single-channel approach has the advantage that no satellite positions need to be known per se and thus no complete navigation messages need to be read and used. Additionally, such an approach also makes the method flexible for processing data from any other (future) GNSS or for parallel processing, which could prove relevant when considering a large number of receivers.
We study the carrier phase-slip and code-outlier detection capabilities of the single-receiver, single-channel model. For the integrity monitoring of carrier-phase data, various studies can already be found in the literature. For dual-frequency GPS data, for instance, methods of carrier phase-slip detection are discussed and tested in (Lipp and Gu 1994;Mertikas and Rizos 1997;Blewitt 1998;Teunissen 1998a;Gao and Li 1999;Jonkman and de Jong 2000;Bisnath and Langley 2000;Bisnath et al. 2001;Liu 2010;Miao et al. 2011). More recent studies on triple-frequency carrier phase-slip detection can be found in (Fan et al. 2006;Dai et al. 2009;Wu et al. 2010;De Lacy et al. 2011;Xu and Kou 2011;Fan et al. 2011).
Our contribution differs from these previous studies, because of its focus on the detectability of single-receiver, singlechannel modeling errors. Next to the phase-slip detection, the detectability of code-outliers is studied as well. Our analysis is analytical, while supported by numerical results. Analytical expressions are derived for the minimal detectable biases (MDBs) of the uniformly most powerful invariant tests (Baarda 1967(Baarda , 1968Teunissen 1990a). The MDB is an important diagnostic tool for inferring the strength of model validation. Examples of such studies for geometry-dependent and integrated GNSS models can be found in (Salzmann 1991;Teunissen 1998b;de Jong and Teunissen 2000;Hewitson and Wang 2010).
This contribution is organized as follows: In Sect. 2 we formulate the multi-frequency, single-receiver, geometry-free GNSS model. This is done for an arbitrary number of frequencies. An overview of the model's redundancy for different measurement scenarios gives a first indication of the model's testability. In Sect. 3 the uniformly most powerful invariant test-statistics for spikes and slips are developed. It is shown how they can be applied to test for code-outliers, phase-slips and ionosphere disturbances. The strength of these test-statistics is described by their corresponding MDBs, for which lower bounds and upper bounds are also given. Due to the relatively simple structure of the geometry-free model, the expression for the MDB can be decomposed into a time-dependent and a time-invariant component. The effect of the time-dependent component is shown in Sect. 3, while the characteristics of the timeindependent part are studied in the sections following. The detectability of phase-slips is treated in Sect. 4. An analytical expression for the phase-slip MDB is derived and it is used to assess the single-, dual-and multi-frequency phase-slip detectability for GPS and Galileo. To evaluate the influence of the code data, the analysis is performed for both the case with code data present and without code data present. The latter case is also of interest, for instance, when one wants to avoid the use of multipath corrupted code data. In Sect. 5, an analytical expression for the codeoutlier MDB is derived. It is used to study the code-outlier detectability for the single-, dual-and multi-frequency GPS and Galileo case, including the case that phase data are absent. In Sect. 6, the MDB for an ionospheric disturbance is presented and analyzed. Finally, the detectability of temporary loss-of-lock on all phase observables is treated in Sect. 7.
2 The multi-frequency, single-receiver geometry-free model

Functional model
The carrier phase and pseudorange (code) observation equations of a single receiver that tracks a single satellite on frequency f j = c/λ j (c is speed of light, λ j is jth wavelength and j = 1, . . . , n) at time instant t (t = 1, . . . , k), are given as (Teunissen and Kleusberg 1998;Misra and Enge 2001;Hofmann-Wellenhoff and Lichtenegger 2001;Leick 2004), where φ j (t) and p j (t) denote the single receiver observed carrier phase and pseudorange, respectively, with corresponding zero mean noise terms n φ j (t) and n p j (t). The unknown parameters are ρ * (t), I(t), b φ j and b p j . The lumped parameter ρ is formed from the receiver-satellite range ρ(t), the receiver and satellite clock errors, cδt r (t) and cδt s (t), respectively, and the tropospheric delay T (t). The parameter I(t) denotes the ionospheric delay expressed in units of range with respect to the first frequency. Thus for the f j -frequency pseudorange observable, its coefficient is given as μ j = f 2 1 / f 2 j . The GPS and Galileo frequencies and wavelengths are given in Table 1. The parameters b φ j and b p j are the phase bias and the instrumental code delay, respectively. The phase bias is the sum of the initial phase, the phase ambiguity and the instrumental phase delay.
Both b φ j and b p j are assumed to be time-invariant. This is allowed for relatively short time spans, in which the instrumental delays remain sufficiently constant (Liu et al. 2004). The time-invariance of b φ j and b p j implies that only time-differences of ρ * (t) and I(t) are estimable. We may therefore just as well formulate the observation equations in time-differenced form. Then the parameters b φ j and b p j get eliminated and we obtain assume information about the absolute ionospheric delays, but rather on the relative, time-differenced, ionospheric delays. We therefore have the additional (pseudo) observation equation with the (pseudo) ionospheric observable I o (t, s). The sample value of I o (t, s) is usually taken to be zero.
and g(t, s) = g(t) − g(s), then the expectation E of the 2n + 1 observation equations of (2) and (3) can be written in the compact vector-matrix form with e n the n-vector of ones and μ = (μ 1 , . . . , μ n ) T . This two-epoch model can be extended to an arbitrary number of epochs. Let y = (y(1) T , . . . , y(k) T ) T and g = (g(1) T , . . . , g(k) T ) T , and let D k be a full rank k × (k − 1) matrix of which the columns span the orthogonal complement of e k = (1, . . . , 1) T , D T k e k = 0. Then dy = (D T k ⊗ I 2n+1 )y and dg = (D T k ⊗ I 2 )g are the time-differenced vectors of the observables and parameters, respectively, and the k-epoch version of (4) can be written as where ⊗ denotes the Kronecker product. The Kronecker product of an a × b matrix M = (m i j ) and a c × d matrix For properties of the Kronecker product, see, e.g., Rao (1973). Model (6), or its two-epoch variant (4), is called the time-differenced, single receiver geometry-free model. It will be referred to as our null hypothesis H 0 .

Stochastic model
The n × n variance matrices of the undifferenced carrier phase and (code) pseudo range observables φ(t) and p(t) are denoted as Q φφ and Q pp , respectively. We assume these variance matrices to be time-invariant and we also assume cross-correlation between phase and code to be absent. Thus for the dispersion of the two-epoch model (4) we have where the scalar σ 2 dI denotes the variance of the timedifferenced ionospheric delay.
To determine the variance matrix of the time-differenced ionospheric delays, let D(I) = Q II be the variance matrix of the absolute ionospheric delay vector I = (I(1), . . . , I(k)) T . The variance matrix of the timedifferenced ionospheric delay vector d I = (D T k ⊗ 1)I is then given as D(dI) = D T k Q II D k . It is through the choice of Q II that we can model the time-smoothness of the ionospheric delays. If we assume that the time series of ionospheric delays can be modeled as a first-order autoregressive stochastic process, then the covariance between I(t) and I(s) is given as σ 2 I β |t−s| , with 0 ≤ β ≤ 1. The two extreme cases are β = 0 and β = 1. In the first case, Q II is a scaled unit matrix and I(t) is considered a white noise process. In the second case, the variance matrix equals the rank-one matrix Q II = σ 2 I e k e T k and I(t) is considered a random constant. In the first case we have D(dI) = σ 2 I D T k D k , while in the second case we have D(dI) = 0.
For the two-epoch case of (7), the variance of the time-differenced first-order autoregressive ionospheric delay works out as For two successive epochs we have σ 2 dI = 2σ 2 I (1 − β), while for larger time-differences the variance will tend to the white-noise value σ 2 dI = 2σ 2 I if β < 1. Thus σ 2 I and β can be used to model the level and smoothness of the noise in the ionospheric delays. We used the above stochastic model for both our analytical and numerical analyses. We have used values for the measurement precision of the multi-frequency GNSS signals reported by Simsky et al. (2008) and de Bakker et al. (2009 which are based on real measurements. The precision of the Galileo E1 signal reported by these publications are in close agreement, for the E5a signal the more conservative value of Simsky et al. (2008) has been adopted. For the GPS L2 signal we will use the same value as for the GPS L1 signal. All zenith-referenced values are summarized in Table 2. To obtain the standard deviations for an arbitrary elevation, these values still need to be multiplied with an elevation dependent function. Several authors studied this dependence, either as function of signal-to-noise (SNR) ratios (e.g., carrier-to-noise density ratio C/N 0 ) or as function of elevation itself, e.g., (Euler and Goad 1991;Ward 1996;Langley 1997;Hartinger and Brunner 1999;Collins and Langley 1999;Wieser et al. 2005). Such weighting will also help suppressing the effect of multipath (de Bakker 2011). For our purposes of studying and evaluating the MDBs, the differences between these functions are negligible. The simplest function, being the cosecant as function of elevation, has a value of about 4 at 15 degrees elevation and reaches its minimum of 1 at 90 degrees elevation.

Redundancy
A prerequisite for being able to perform statistical tests is the existence of redundancy. For a full rank model, redundancy is defined as the number of observations minus the number of unknown parameters. We have summarized the redundancy of the k-epoch model (6) in Table 3. We discriminate between the ionosphere-weighted case and the ionosphere-float case. In the latter case, no a priori information is assumed about the ionospheric delays. Hence, in this case, all ionospheric delays are treated as completely unknown. This results therefore in a redundancy reduction of k − 1, being the number of unknown time-differenced ionospheric delays. We also discriminate between the phase and code case, and the code-only (phaseless) and phase-only (codeless) cases. When both phase and code data are used, the ionosphere-weighted redundancy equals (k − 1)(2n − 1). Thus in this case, redundancy exists for any number of frequencies, provided k ≥ 2. That at least two epochs of data are needed is of course due to the fact that we are working with time-differenced data. That already single-frequency (n = 1) processing provides redundancy is due to the ionospheric information. Without this information, there would be no redundancy in the single-frequency case, but only in the dualand multi-frequency cases, provided both phase and code data are used. The phase-only and code-only redundancies are the same. In the phase-only and code-only cases we have (k − 1)n observations less than in the phase and code case. Hence, this is the number by which the redundancy drops when either the code data or the phase data are left out. Thus in the phase-only or code-only cases, single-frequency testing is impossible even if ionospheric information is provided.

Testing and reliability
In this section we formulate our alternative hypotheses and present the corresponding test statistics.

Outliers, cycle slips and loss-of-lock
We now formulate our alternative hypotheses for the singlereceiver, geometry-free GNSS model. They accommodate model biases such as outliers in the pseudo range data, slips in the carrier phase data and loss-of-lock.
Recall that the undifferenced observational vector of epoch t is given as y(t) = (φ(t) T , p(t) T , I(t)) T . Now assume that a model error has occurred in the data of epoch l and that this (2n + 1)-bias vector can be parametrized as Hb, where H is a given matrix of order (2n + 1) × q and b is an unknown vector having q entries. Then Hb is the difference between the expectation of y(l) under the null hypothesis H 0 and the expectation of y(l) under the alternative hypothesis H a . Thus E(y(l)|H a ) = E(y(l)|H 0 ) + Hb. Through the choice of matrix H , we can describe the type of model error. For instance, if all the phase data are assumed erroneous, as would be the case after a temporary loss-of-lock on all phase observables, then H = (I n , 0, 0) T and q = n. But if only the pseudo-range data on frequency j is corrupted with an outlier, then H = (0, δ T j , 0) T and q = 1, where δ j is an n-vector having a 1 as its jth entry and zeros elsewhere.
Apart from describing the model error through matrix H , we also need to specify the time behavior of the model error.
Here we consider spikes and slips. A model error behaves as a spike if it occurs at one and only one epoch. A model error is said to behave as a slip if it persists after occurrence. Examples of spikes are outliers in the pseudo-range data or in the ionospheric delays. Examples of slips are cycle slips in the phase data or momentary loss-of-lock.
If we assume the model error Hb to behave as a spike at epoch l, then E(y|H a ) = E(y|H 0 ) + (s l ⊗ H )b, where s l is a k-vector having a 1 as its lth entry and zeros elsewhere. Would we assume the error to be persistent, however, as would be the case after a loss-of-lock or after a slip, then s l is a k-vector having zeros in its first l −1 entries, but 1s in all its remaining entries.
Thus with suitable choices for the vector s l and the matrix H , one can model outliers in the code data, cycle slips in the phase data, disturbances in the ionosphere and even a complete loss-of-lock. The formulation of the alternative hypotheses in terms of the time-differenced data follows then from pre-multiplying The null-and alternative hypotheses treated in the present contribution are therefore given as In order to test H 0 against H a , the uniformly most powerful invariant (UMPI) test is used, see e.g., (Arnold 1981;Koch 1999;Teunissen 2000). It rejects the null hypothesis in favor of the alternative hypothesis, if whereb, with variance matrix Qbb, is the least-squares esti- . Hence, with α being the probability of wrongful rejection, the critical value χ 2 α (q, 0) of test (12) is

Test statistics for spikes and slips
In order to derive the appropriate test statistics, we first determine the least-squares estimator of b in (9). Here and in the following we assume D(dI) = σ 2 I D T k D k and therefore D(y(i)) = blockdiag(Q φ , Q p , σ 2 I ) call = Q. The leastsquares estimator of b and its variance matrix is given in the following theorem: , the least-squares estimator of b under H a (cf. 9) and its variance matrix are given aŝ and the uniformly most powerful invariant test statistic for testing H 0 against H a is given as The bias estimator and its variance matrix are given the arguments l and k to emphasize that it is an estimator of an error occurring at epoch l, based on k epochs of data.
From the Kronecker product structure of (13) it follows thatb(l, k) and its variance matrix can be computed directly from its single-epoch counterparts. We therefore have the following result:

Corollary 1 Let the single-epoch bias estimator and its variance matrix be given aŝ
and define the (2n and Note that the time dependency in the coefficients of (17) is captured bys l , which will be different for spikes and slips.
Spikes For spike-like biases the k-vector s l is a unit vector having a 1 as its lth entry. If we make use of This last expression clearly shows thatb(l, k) is the difference of the local estimatorβ(l) and the time-average of theβ(i) over the error-free time instances i = 1, . . . , k; i = l. Under H a the mean of the local estimator is b and that of the timeaverage is zero. Although one can use the local estimator β(l) as the bias-estimator, the local estimator does not make use of the information that all other epochs are assumed to be bias-free. Hence, the local estimator can be improved by including this information through subtraction of the zeromean time-average of the bias-free time instances.
For computational purposes it is easier to use the first expression of (18), because of the presence of the running average (which can also be computed recursively). We therefore use this expression to formulate the corresponding test statistic for spikes. It reads with the running averageβ(k) = 1 The test statistic (19) can be used for any l ≤ k. However, when l < k, there is a delay in testing of k − l. For some applications this may not be acceptable or necessary. For those cases one will use (19) with l = k, which gives Thus in this case, the test statistic is formed from the difference of the local bias estimatorβ(i) for i = k and its time-average over the previous k − 1 epochs,β(k − 1). Slips For slip-like biases the k-vector s l is a vector having 1s as its last k − l + 1 entries and zeros elsewhere. If we make use of This last expression clearly shows that in case of a slip, b(l, k) is the difference of two time-averages of theβ(i). The first time-average averages over all epochs that supposedly contains the bias, while the second is the time-average over all bias-free epochs. Under H a , the mean of the first time-average is b, while that of the second time-average is zero.
Note that (21) reduces to (18) when l = k. This shows that one cannot discriminate between spikes and slips on the basis of one single epoch. That is, one needs to have a delay (k > l) to be able to separate spikes from slips.
For computational purposes it is easier to use the first expression of (21), because of the presence of the running averages. We therefore use this expression to formulate the corresponding test statistic for slips. It reads It is thus formed from the difference of two time-averages of theβ(i), namely the time-average over the epochs up to and including the time of testing k and the time-average over the epochs up to time l at which the slip is assumed to have started. For l = k this test statistic becomes identical to (20).

Minimal detectable biases
Under the alternative hypothesis H a , the test statistic T q is distributed as a noncentral χ 2 -distribution with q degrees of (12) is an UMPI-test, meaning that for all b, it maximizes the power within the class of invariant tests. Here, power, denoted as γ , is defined as The power of test (12) depends on the degrees of freedom q (i.e., the dimension of b), the level of significance α, and through the noncentrality parameter λ, on the bias vector b. Once q, α and b are given, the power can be computed.
One can, however, also follow the inverse route. That is, given the power γ , the level of significance α and the dimension q, the noncentrality parameter can be computed, symbolically denoted as λ 0 = λ(α, q, γ ). Figure 1 shows the relation between λ 0 and γ for different values of q and α. With λ 0 given, one can formulate the following quadratic This equation is said to describe, for test (12), the reliability of the null hypothesis H 0 with respect to the alternative hypothesis H a . For q = 1, Eq. (23) describes an interval, for q = 2 it describes an ellipse and for q > 2 it describes an (hyper)ellipsoid. Bias vectors b ∈ R q that lie on or outside the ellipsoid (23) can be found with at least probability γ .
To determine the bias vectors that satisfy (23), we use the factorization b = ||b||d, where d is a unit vector (d T d = 1). Substitution into (23) and inversion gives This is the celebrated Minimal Detectable Bias (MDB) vector of Baarda's reliability theory (Baarda 1967(Baarda , 1968). Its length is the smallest size of bias vector that can be found with probability γ in the direction d with test (12). By letting d vary over the unit sphere in R q one obtains the whole range of MDBs that can be detected with probability γ with test (12). For the one-dimensional case (q = 1), we have d = ±1 and therefore MDB(l, k) = σb(l, k) √ λ 0 . Baarda, in his work on the strength analysis of general purpose networks, applied his general MDB-form to data snooping, thus obtaining the scalar boundary values (in Dutch: 'grenswaarden'). Applications of the vectorial form can be found, for example, in Van Mierlo (1980 and Kok (1982a) for deformation analysis, in Foerstner (1983) for photogrammetric linear trend testing and in Teunissen (1986) for testing digitized maps. For recursive testing, with application of testing time and types of error, the first innovation-based vectorial MDB was given in Teunissen (1990b). Application of the generalized eigenvalue problem to the vectorial form to obtain MDB-bounds can be found in e.g., (Teunissen 2000;Knight et al. 2010).
Earlier (cf. 12) it was assumed that the variance matrix Qbb in the teststatistic T q is known. In case, however, the variance factor σ 2 in Qbb = σ 2 Gbb is unknown, then the Chisquare distributed teststatistic needs to be replaced by the F-distributed teststatistic T q,d f =b T G −1 bbb /(σ 2 q), in whicĥ σ 2 is the unbiased estimator of the variance factor under the alternative hypothesis and d f is its degrees of freedom. This teststatistic is distributed under H a as T q,d f ∼ F(q, d f, λ), with the same noncentrality parameter as that of T q (Kok 1982b;Koch 1999). Therefore also in case σ 2 is estimated, will the same MDB expression be found as given in (24). Hence, similarly as with the other statistical parameters, the effect on the MDB, for the two cases σ 2 known versus σ 2 estimated, is only felt through the scaling factor λ 0 . In this contribution we work with λ 0 computed from the Chi-square distribution.
If we substitute (16) into (24), the MDBs for spikes and slips work out as and The MDBs get smaller if more epochs of data are used (k gets larger) and/or a more precise bias-estimatorβ is used (Qββ gets 'smaller'). Note that the spike-MDB depends on k, whereas the slip-MDB depends on both l and k. Hence, the detectability of spikes does not depend on the time instance  (k = 4, 5, 6, 8, 10, 20) of their occurrence, whereas the detectability of slips does depend on these time instances. For k fixed, the slip-MDB reaches its minimum at k/2 + 1. Hence, for a given observational time span, slips are best detectable if they occur 'half way' in the time span. The graph of the function f (l, k) is given in Fig. 2. It shows by how much the MDB(2, 2) improves when l and k are larger than 2.
In the following, we present the MDBs for the two-epoch case (k = 2, l = 2). The corresponding MDB-values for arbitrary epochs can then be obtained by using the multiplication factor f (l, k) of (25).

Detectability of carrier phase slips
In this section we present and analyze the phase-slip MDBs. This is done for the single-, dual-and multi-frequency case. We also analyze the phase-slip MDBs in case code data are absent.

Minimal detectable carrier phase slips
In the following, MDB values have been computed numerically using the complete set of zenith-referenced code and phase variances according to Table 2, and analytical MDB expressions have been derived for which the phase and code variance matrices have been simplified to scaled unitmatrices, i.e., Q φφ = σ 2 φ I n and Q pp = σ 2 p I n . For the standard deviation of the scaled unit-matrices we have used the mean value of the standard deviations of the available signals (except when explicitly stated otherwise). The MDB values from both the numerical computations and the analytical expressions are presented graphically. The differences between the two are shown to be small, especially when the used signals have comparable precision, which indicates that the analytical expressions indeed give a proper representation of the single-receiver, single-channel detectability (in the MDB graphs below, the dashed curves correspond to the analytical expressions, while the full curves correspond to the MDBs computed from the actual variance matrices). For all MDBs which pertain to an error on a single frequency, we have chosen to compute the MDB for an error on the first frequency of the available frequencies. For GPS this is the L1 frequency and for Galileo the E1 frequency when available.
The analytical expression for the multi-frequency phaseslip MDB is given in the following theorem: φ I n and Q pp = σ 2 p I n . Then for k = 2, the MDB for a slip in the jth frequency carrier phase observable is given as with scalar where = σ 2 φ /σ 2 p is the phase-code variance ratio andμ = 1 n n i=1 μ i is the average of the squared frequency ratios Proof see Appendix.
Note that the slip MDB is scaled with σ φ , the small standard deviation of the carrier phase observable. Thus if the bracketed ratio of (28) is not too large, small carrier phase slips will be detectable. The bracketed term depends on λ 0 and on n * . The scalar n * is dependent on , μ i (i = 1, . . . , n) and σ 2 φ /σ 2 dI . Hence, it depends on the measurement precision, on the number of frequencies and their spacings, and on the smoothness of the ionosphere. The MDB gets smaller if gets larger, i.e., if more precise code data are used. The MDB also gets smaller if n gets larger, i.e., if more frequencies are used. Finally, note that the MDB-dependence on the frequency diversity (i.e., on μ i , i = 1, . . . , n) is driven to a large part by the smoothness of the ionosphere. This dependence is absent for σ dI = 0 and it gets more pronounced the larger σ dI gets.
We now analyze the slip MDB for the GPS and Galileo single-, dual-and triple-frequency cases.

Single-frequency case
For the single-frequency case (n = 1), the carrier phase slip MDB (28) can be worked out to give This result shows that single-frequency phase slip detection is possible in principle, but that its performance depends very much on the smoothness of the ionosphere and on the measurement precision of code. If σ 2 dI = ∞, then MDB φ j = ∞, meaning that single-frequency slip detection has become impossible. If the ionospheric delay is so smooth that the ratio σ 2 dI /σ 2 p may be neglected, we get with ≈ 0, MDB φ 1 = σ p √ 2λ 0 . Hence, for a sufficiently smooth ionospheric delay, we may detect slips of about six times the code standard deviation (here and in the following the reference value of the noncentrality parameter is taken as λ 0 = 17.02; it corresponds to q = 1, α = 0.001 and γ = 0.80, see Fig. 1).
The numerical values of the GPS and Galileo singlefrequency phase-slip MDBs are graphically displayed in Fig. 3 as function of σ dI . It shows that the MDBs are approximately constant for σ dI ≤ 10 −2 m, but then rapidly increase for larger values of σ dI . The constant values are about 146 cm for L1 and L2, 117 cm for E1, 88 cm for L5, E5a, E5b and E6 and 41 cm for E5. These three levels reflect the difference in the code measurement precision of the signals. Since these single-frequency MDBs are all larger than their corresponding wavelengths, time-windowing (N = k − l + 1 > 1 c.f. 25) will be needed to bring the MDBs down to smaller values. Accepting a delay in testing is then the price one has to pay for the increase in detectability.

Dual-and multi-frequency case
It follows from (28) that the dual-and multi-frequency phaseslip MDB is bounded from below as This lower bound is obtained for σ dI = 0. This shows that for a sufficiently smooth ionospheric delay, very small slips can be found, even for the two-epoch case. This is confirmed by the S-shaped MDB graphs of Fig. 4. Very small slips can be detected (in order of a few cm), if the ionospheric delays are sufficiently smooth (σ dI ≤ 10 −2 m). In this case, there is also no significant difference between the dual-frequency and multi-frequency performance. This difference is present though, for larger values of σ dI . In particular the dual-frequency GPS phase-slip MDBs, and the Galileo E1 E5a, increase steeply when σ dI ≥ 10 −2 m gets larger. For the multi-frequency phase-slip MDBs the increase is much more moderate. All the multi-frequency phase-slip MDBs remain smaller than 7 cm, whereas the dual-frequency MDBs are all smaller than 22 cm. Those of Galileo are smaller than their GPS counterparts, due to their improved code precision, except that L1-L2-L5 outperforms E1-E5a-E5b due to a better distribution of the frequencies. We can also see that the analytical expression is further removed from the numerical results for the E1-E5 combination, which is a result of the large difference in precision between these two signals. The scaled unit matrix used for the code variance approximates the actual variance matrix less well.  Since code data are generally much less precise than carrier phase data, one may wonder whether code data are really needed for carrier phase slip detection. This is also of interest for those applications where the code data are corrupted by multipath. In this section we therefore investigate what happens when σ p → ∞. The corresponding MDBs can be obtained from Theorem 1 by taking the limit lim σ p →∞ MDB φ j . We have the following result.
Lemma 1 (Codeless phase-slip MDB) The codeless phaseslip MDB lim σ p →∞ MDB φ j follows from (28) using or using the limit of the inverse whereμ ( j) = 1 n−1 n i = j μ i is the average of the squared frequency ratios that excludes μ j .
The above two expressions clearly show the effect of frequency diversity. The MDB gets smaller, the larger the frequency diversity n i=1 (μ i −μ) 2 is. And within a set of n > 1 MDBs, the smallest MDB φ j is the one for which μ j is closest to the averageμ.

Single-frequency case
In the single-frequency case we have n = 1 and thus no frequency diversity, i.e., n i=1 (μ i −μ) 2 = 0. This shows that n * = 1 for n = 1 (c.f. 32) and that therefore MDB φ j = ∞. Hence, phase-slip detection without code data is impossible for the single-frequency case (see also the redundancy Table 3).

Dual-frequency case
In the dual-frequency case (n = 2) we have (μ j −μ) 2 = 1 2 ( n i=1 (μ i −μ) 2 ). Hence, n * = 1 (c.f. 32) if n = 2 and σ dI = ∞. In that case dual-frequency phase-slip detection without code data becomes impossible. In all other cases, however, we have This shows that very small phase-slips can be detected when the ionospheric delays are sufficiently smooth. Figure 5 shows them to be less than 3 cm for σ dI ≤ 10 −2 m. For larger values, the MDBs increase linearly, with the gradient driven by the frequency diversity; the closer the frequencies, the less steep the increase is. When we compare Fig. 5 with Fig. 4, we note that the phase-slip MDB values are not too different for sufficiently small σ dI , but that their differences increase for larger σ dI . Thus the presence of the code data are particularly needed when the ionospheric delays are not smooth enough. Codeless dual-frequency phase-data is sufficient to detect phaseslips otherwise.

Multi-frequency case
In the multi-frequency case the general form of the codeless phase-slip MDB follows from (28) and (33) as To discuss its dependence on the ionospheric variance and on the frequency distribution, we start from the most optimal scenario. The MDB is smallest when σ 2 dI = 0, in which case the frequency dependence only includes the number of frequencies but not their diversity, For σ 2 dI = 0, the MDB is smallest if all frequencies are the same (μ i = 1, i = 1, . . . , n), i.e., if the vector μ = (μ 1 , . . . , μ n ) T is parallel to the vector e n = (1, . . . , 1) T . One would then get the same MDB as (36), i.e., the one that corresponds with σ 2 dI = 0. This can be explained as follows: If μ and e n are parallel vectors, then the ionospheric delay I(i) gets lumped with ρ * (i) (c.f. 4 and 5), reducing the model in essence to an ionosphere-fixed one. Thus in the absence of code data, the absence of frequency diversity is best for phase-slip detection.
Now assume that not all components of vector μ are the same, but that instead all but one of them are the same. Thus μ i =μ, ∀i = j. Then the codeless phase-slip MDB (35) works out as This expression generalizes (34) for n > 2. Its value goes to infinity for σ 2 dI → ∞. In this case, it is the lack of frequency diversity in the n − 1 phase data, φ i , i = j, that makes it impossible to solve for the ionospheric delay and therefore also for detecting a slip in the jth frequency phase observable.
For the triple-frequency case (n = 3), the following bounds can be formulated: with the cyclic indices ( j 1 , j 2 , j 3 ) = (1, 2, 3), (2, 3, 1), (3, 1, 2) depending on whether j = 1, 2 or 3. The lower bound corresponds with σ 2 dI = 0, while the upper bound corresponds with σ 2 dI = ∞. These bounds show that very small phase-slips can be detected, even when the code data are absent. This is confirmed by the graphs of Fig. 6. All codeless phase-slip MDBs are less than 10 cm and even as small as about 1 cm when σ dI ≤ 3 × 10 −3 m. For larger values of σ dI the differences in the triple-frequency MDBs become clearly visible. This is due to their frequency dependence.
When we compare the codeless results of Fig. 6 with the triple-and quadruple-frequency results of Fig. 4, no big differences can be seen. Hence, the impact of the code data on the phase-slip MDBs becomes less pronounced if more than two frequencies are used. The only noteworthy difference between the results of the two figures is the performance of E1-E5a-E5b. This difference in performance, compared with the other triple-frequency results, is due to the small frequency separation of E5a and E5b.

Detectability of code outliers
In this section we present and analyze the code-outlier MDBs. This is done for the single-, dual-and multi-frequency case. We also analyze the code-outlier MDBs in case phase data are absent.

Minimal detectable code outliers
The analytical expression for the multi-frequency codeoutlier MDB is given in the following theorem: Theorem 3 (Code-outlier MDB) Let H = (0, δ T j , 0) T , Q φφ = σ 2 φ I n and Q pp = σ 2 p I n . Then for k = 2, the MDB for an outlier in the jth frequency code observable is given as where = σ 2 φ /σ 2 p is the phase-code variance ratio andμ = 1 n n i=1 μ i is the average of the squared frequency ratios Proof As the phase observables and code observables play a dual role in the two-epoch model (4), the code-outlier MDB can be found from the expression of the phase-slip MDB (28) through an interchange of the phase and code variance.
Although (39) and (40) have the same structure as (28) and (29), respectively, there are marked differences between these expressions. First note that (39) is scaled by the standard deviation of code and not by that of phase as in (28). Second, the frequency-dependent term between brackets in (40) is multiplied with the very small phase-code variance ratio , while this is not the case with the bracketed term of (29). The consequences of these differences are that the dual-and multi-frequency code-outlier MDBs are generally larger than those of the phase-slip MDBs and that they are less sensitive to the frequencies. The exception occurs in the single-frequency case.

Single-frequency case
In the single-frequency case, for k = 2, the code-outlier MDB is identical to that of the phase-slip MDB. This follows if we interchange the role of σ 2 p and σ 2 φ in (30). For k > 2, however, the MDBs differ of course. In case of a slip the multiplication factor is 1 √ 2 1 k−l+1 + 1 l−1 , while for an codeoutlier it is 1 √ 2 1 + 1 k−1 .

Dual-and multi-frequency case
We already remarked that in case of a sufficiently smooth ionosphere, the dual-and multi-frequency phase-slip MDBs become relatively insensitive to the frequencies. This can be seen in Fig. 4, but it also follows from the lower bound (31). The reason for this insensitivity lies in the fact that the frequency-dependent bracketed term of (29) reduces to 1 for σ dI → 0. This effect is also present in the code-outlier MDB bracketed term of (40). However, with (40) there is the additional effect that the bracketed term is multiplied with the very small phase-code variance ratio. Hence, the MDB is now also less sensitive to the frequencies for larger values of σ dI . This means that one can expect the lower bound to be a good approximation to the code-outlier MDB. After all, this lower bound follows from neglecting 1 m * in (40). Thus, for α = 0.001 and γ = 0.80, giving λ 0 = 17.02, one can expect the code-outlier MDB to be about six times the code standard deviation. This is confirmed by the graphs of Fig. 7. The figure shows that the code-outlier MDBs are not very sensitive to the available signals and frequencies. Additionally, the MDBs are nearly constant with the variance of the time-differenced ionospheric delay (no MDB variations are visible at the scale of this figure). Compare with phaseslip MDB Fig. 4. The two levels shown in Fig. 7 correspond to the different code measurement precision levels of the GPS L1 and Galileo E1 signals (cf. Table 2). Figure 7 is the only figure for which we have directly used the standard devia-  Fig. 7 Dual-and multi-frequency GPS and Galileo code outlier MDBs for k = 2; dashed curves are based on Eqs. (39) and (40) tion of the L1 and E1 signals in the analytical expressions (instead of the mean value of the available signals), since only the precision of these signals has any impact on the size of the MDBs.

Code outlier detection without phase data
The lower bound (41) follows from taking the limit σ 2 φ → 0 in (39). We now consider the other extreme, σ 2 φ → ∞. It corresponds to code outlier detection without the use of carrier phase data. The corresponding MDB will be larger than (41). Furthermore, the absence of carrier phase data will now make the MDB dependent on the frequencies and the ionospheric variance. This follows, since the phase-variance limit of (40) reads We now again discriminate between the three cases n = 1, n = 2 and n > 2.

Single-frequency case
Just like single-frequency codeless phase-slip detection is impossible, so is single-frequency code-outlier detection without phase data. This follows directly from Table 1, which shows that redundancy is absent, when phase data is absent in case n = 1. It also follows from (42), which shows that m * = 1 if n = 1.

Dual-frequency case
In the dual-frequency case, the phaseless code-outlier MDB is given as It follows from interchanging the role of the phase-and code variance in (34). Expression (43) shows that the MDB gets larger (poorer detectability) if the frequency separation gets larger (larger |μ 2 −μ 1 |). This is perhaps contrary to what one would expect. However, one should not confuse the estimability of the ionosphere in the absence of biases, i.e., under the null hypothesis H 0 , with the detectability of biases in the presence of the ionosphere. Since the variance of the ionospheric delay can be shown to be inversely proportional to |μ 2 −μ 1 |, the ionospheric delay estimator gets indeed more precise when the frequencies are further separated. The contrary happens, however, with the outlier detectability. The detectability of the outliers becomes poorer for larger frequency separation and this effect becomes more pronounced for larger σ 2 dI . This frequency dependence is absent in case σ 2 dI = 0. The graphs for the dual-frequency phaseless code-outlier MDBs are given in Fig. 8. Compare this figure with Fig. 5. The graphs in both figures show a similar shape. In the case of the code-outlier MDBs, however, the graphs do not coincide because of the different levels of code measurement precision of the signals. When we compare Fig. 8 with Fig. 7, we also clearly show the impact of the phase data. With the phase data included, the code-outlier MDB not only becomes smaller, but also the steep increase experienced in the phaseless case for σ dI ≥ 10 −1 m is eliminated. Without the phase data, the  39) and (42) code-outlier MDB is more dependent on the smoothness of the ionospheric delay.

Multi-frequency case
The multi-frequency phaseless code-outlier MDB follows from interchanging the phase and code variance in (35). Their graphs are shown in Fig. 9. When compared with Fig. 8, they seem to show the same behavior as in the dual-frequency case. This is not true, however. In the dual-frequency case the code-outlier MDB keeps growing with increasing σ dI , whereas in the multi-frequency case the MDB levels off for large enough σ dI . Hence, for larger values than σ dI = 10 0 m, the graphs of Fig. 9 are S-shaped, just like those of Fig. 6.

Detectability of ionospheric disturbances
The analytical expression for the multi-frequency ionospheric disturbance MDB is given in the following theorem: Theorem 4 (Ionospheric MDB) Let H = (0, 0, 1) T , Q φφ = σ 2 φ I n and Q pp = σ 2 p I n . Then for k = 2, the MDB for an ionospheric disturbance is given as where = σ 2 φ /σ 2 p is the phase-code variance ratio,μ = 'variance' of the squared frequency ratios Proof LetdI, with variance (45), be the LS estimator of dI based on the two-epoch model under H 0 (cf. 9) using phase and code data only. Then it follows from the structure of the model under H a that the LS bias estimator of the ionospheric disturbance is given as the differenceb dI = dI −dI. Therefore, its variance is given by the sum σ 2 dI + σ 2d I , from which the result follows.
The above result clearly shows what role is played by the a-priori ionospheric standard deviation, σ dI , the measurement precision, σ φ and σ p , and the distribution of the frequencies, f i , i = 1, . . . , n. In case all n frequencies are equal, then σ 2 μ = 0, and the MDB reduces to Although ionospheric disturbance testing is then still possible in principle, its performance will then primarily be driven by the code variance. Figure 10 shows the MDB for the singlefrequency case (n = 1). When there is a nonzero 'variance' in the frequencies, σ 2 μ = 0, then codeless or phaseless testing is possible as well. The codeless MDB follows from (44) as If we replace σ 2 φ in this expression by σ 2 p , we obtain the corresponding phaseless MDB. Figures 11 and 12 show the codeless and phaseless MDBs for the different cases.

Detectability of phase loss-of-lock
Phase loss-of-lock is defined as the simultaneous occurrence of an unknown multivariate slip in all n carrier phase observables. To study its detectability, we first determine the variance matrix of the multivariate slip under H a and then provide bounds on the norm of the phase loss-of-lock MDBvector.

The variance matrix of the multivariate slip
In the presence of a phase loss-of-lock, the design matrix of the alternative hypothesis (9) takes for k = l = 2 the form [G, H ], with H = [I n , 0, 0] T . From the structure of [G, H ], it follows that the carrier phase vector φ(t, s) will not contribute to the estimation of the parameters ρ * (t, s) and I(t, s) under H a . These parameters are therefore solely determined by the code observables and a priori ionospheric information. As a consequence, the two-epoch bias estimator is given as the differenceb = φ(t, s) −φ(t, s), whereφ(t, s) = e nρ * (t, s) − μÎ(t, s) is the least-squares phase estimator based solely on the code observables and a priori ionospheric information. Solving forρ * (t, s) and I(t, s), followed by applying the variance propagation law tob = φ(t, s) − e nρ * (t, s) + μÎ(t, s) gives then the variance matrix of the multivariate slip. The result is given in the following Lemma: Lemma 2 (Variance matrix of multivariate slip) For k = l = 2 and H = (I n , 0, 0) T , the variance matrix of the leastsquares estimator of b in (9) is given as with σ 2 dÎ = ( 1 2 μ T Q −1 pp P ⊥ e n μ+σ −2 dI ) −1 , P ⊥ e n = I n − P e n , P e n = e n (e T n Q −1 pp e n ) −1 e T n Q −1 pp and R e n = I n + P e n .
Note that the variance matrix is a sum of three matrices, the entries of which may differ substantially in size. The first matrix is governed by the precision of the phase observables and will therefore have small entries. The second matrix is governed by the precision of the code observables and will therefore have generally much larger entries than the first matrix in the sum. The third matrix depends, next to the precision of the code observables, also on μ and σ 2 dI . Its entries will become smaller if σ 2 dÎ gets smaller. This happens for smaller σ 2 dI (smoother ionospheric delays) and/or larger μ T Q −1 pp P ⊥ e n μ (better code precision and/or larger frequency diversity). Thus if σ 2 dI = ∞, frequency diversity is needed (i.e., μ T Q −1 pp P ⊥ e n μ = 0) so as to avoid the entries of the third matrix in (48) to become infinite.
We now discuss the detectability of the phase loss-of-lock for n ≥ 2. The single-frequency case n = 1 is already treated (cf. 30), since phase loss-of-lock is then equivalent to having a phase-slip.

Bounding the MDB vector
The variance matrix (48) can be used together with (25) to determine the phase loss-of-lock MDB-vector. Its length depends on the direction d in which the multivariate slip occurred, For certain directions, this may result in a short vector, while for other directions, it may be a long vector.
For its length, the following upper bound can be used: where σ 2 d Tb = d T Qbbd is the variance of d Tb , with d being the unit vector that points in the direction of the slip. The upper bound is easier to obtain as it avoids the inversion of the variance matrix (48).
Considering (48), it is not difficult to see in which directions the MDB-vector will be short. If d ⊥ e n , then P T e n d = 0, meaning that the second term of (48) will not contribute. And if d ⊥ span{e n , μ}, then R T e n d = 0 and P T e n d = 0, meaning that now also the third term will not contribute. Thus This shows that the length of the MDB vector is very small indeed if d lies in the (n − 2)-dimensional space span{e n , μ} ⊥ . In this case, the phase loss-of-lock has very good detectability, since the MDB is then solely driven by the very precise carrier phase data. The situation changes drastically, however, if d is taken to lie in span{e n , μ}. In that case the large code variance dependent second and third term of (48) will contribute as well. The largest possible value that the length of the MDB vector can take corresponds with d being the eigenvector of Qbb with largest eigenvalue. The corresponding bounds for the ionosphere-fixed and ionosphere-float case are given in the following Lemma: Lemma 3 (Phase loss-of-lock MDB upper bounds) Let Q φφ = σ 2 φ I n , Q pp = σ 2 p I n . Then Proof see Appendix.
The ionosphere-fixed upper bound corresponds with a slip in the e n -direction, while the ionosphere-float upper bound corresponds with a slip in the direction e n + √ n ||μ|| μ. Thus phase loss-of-lock is most difficult to detect if it results in a slip with such direction. For example, for the ionospherefixed, dual-and triple-frequency (q = 2, 3) cases, the length of the MDB-vector will then be about six to seven times the code standard deviation.
To conclude, we note that we can apply the duality of phase and code to the above lemma and so obtain directly also the length of the MDB vector ||MDB|| p that corresponds to the case that the complete n-dimensional code vector is outlying.
By interchanging the phase-and code variance in (52), we obtain The ionosphere-fixed upper bound is the same as in (52), but the ionosphere-float upper bound differs. In this second bound we note that the effect of frequency diversity (i.e., the effect of f ) gets reduced due to its multiplication with the small phase-code variance ratio . This corresponds to our earlier findings in Sect. 5.1.2, where it was stated that over the considered range of σ 2 dI , the code-outlier MDB is much more insensitive to frequency diversity than the phase-slip MDB is.

Summary and conclusions
In this contribution we presented an analytical and numerical study of the integrity of the multi-frequency single-receiver, single-channel GNSS model. The UMPI test statistics for spikes and slips are derived and their detection capabilities are described by means of the concept of minimal detectable biases (MDBs). Analytical closed-form expressions of the phase-slip and code outlier MDBs have been given, thus providing insight into the various factors that contribute to the detection capabilities of the various test statistics. This was also done for the phaseless and codeless cases, as well as for the case of loss-of-lock.
The MDBs were evaluated numerically for the several GPS and Galileo frequencies. From these analyses it can be concluded that detectability of dual-and triple-frequency phase-slips works well for k = l = 2. Single-frequency phase-slip detectability, however, is problematic, thus requiring more epochs of data.
From the codeless phase-slip MDBs it follows that detectability is not possible in the single-frequency case, but that it is possible for the dual-and triple-frequency cases. In the dual-frequency case the codeless phase-slip MDBs are less then 10 cm if σ dI ≤ 3 cm, while for the triple-frequency case this is already true for σ dI ≤ 1 m. In the triple-frequency case, the phase-slip MDBs get even as small as a few centimeters if σ dI ≤ 3 cm. These codeless results are important as it shows that in the presence of code multipath, one can do away with the code data and then still have integrity for phase slips.
The code outlier MDBs are, except for the singlefrequency case, all relatively insensitive to the smoothness of the ionosphere. The effect of the frequencies is also hardly present in the multi-frequency code-outlier MDBs. Their value is predominantly determined by the precision of the code measurements. From the phaseless code-outlier MDBs it follows that, except for the single-frequency case, code outlier detection is still possible. The multi-frequency phaseless code-outlier MDBs all lie around 1 m for σ dI ≤ 10 cm. They increase rapidly, however, for larger values of σ dI . Acknowledgments The first author is the recipient of an Australian Research Council (ARC) Federation Fellowship (project number FF0883188). Part of this work was done in the framework of the project 'New Carrier-Phase Processing Strategies for Next Generation GNSS Positioning' of the Cooperative Research Centre for Spatial Information (CRC-SI2). The work of the second author was partly supported by the EU Marie Curie program in the framework of the FP7-PEOPLE-IAPP-2008 (SIGMA) project. All this support is gratefully acknowledged.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.

Proof of Theorem 3 (Phase loss-of-lock MDB upper bounds):
We give the proof for the ionosphere-fixed case. For the ionosphere-float case it goes along similar lines.
To determine the eigenvalues of Qbb (c.f. 48) for Q φφ = σ 2 φ I n , Q pp = σ 2 p I n and σ 2 dI = 0, we need to determine the root γ of |2σ 2 φ I n + 2σ 2 p e(e T e) −1 e T − γ I n | = 0 Let D be an n × (n − 1) basis matrix of the orthogonal complement of e. Then the n × n matrix M = [D, e] is of full rank. Pre-and post-multiplication of the matrix in (58) with M T and M allows to reduce the determinantal equation into a product, This shows that there are (n − 1) eigenvalues with the value γ = 2σ 2 φ and one eigenvalue, the largest, with the value γ = 2(σ 2 φ + σ 2 p ) = 2σ p (1 + ). Hence, for the MDB upper bound we get This concludes the proof.