Practical multivariate statistical multipath detection methods

Fig. 1 Simulated GPS three-frequency (red: L1, blue: L2, green: L5) multipath error in PRN02 in the LCPC dataset (relative permittivity =3.9)

GPS modernization is being undertaken. The GPS Block IIR-M satellites are already transmitting an unencrypted civil signal (L2C) on L2 frequency and Block III satellites will transmit a new civil signal (L1C) on L1. Moreover, the signal power of L2 will be increased. This will make tracking of L2 much easier and more reliable and will increase the use of L2 in high precision kinematic applications. An additional signal, the so-called L5, will be available on GPS Block IIF satellites scheduled for launch beginning in late 2009. Both the modernised L2 and the new L5 civil signals allow coherent tracking of code and phase and so avoid the losses that occur when tracking the current P(Y) code in L2. This had led to the extensive current interest, e.g. (Hatch et al., 2000) in investigating the potential of three-frequency data for a wide range of applications. On the other hand, the European GNSS, named Galileo, is being developed to provide four carrier frequencies and its Full Operational Capability (FOC) is scheduled to be in 2013. Galileo signals are expected to be available to users in four categories: Open Service (OS), Safety-of-Life (SoL) service, Commercial Service (CS), and Public Regulated Service (PRS).

This paper investigates the application of a sophisticated multiple outlier detection technique, which the author refers to as the cocktail multiple outlier detection algorithm, to multipath contaminated measurements in the three-frequency GPS and Galileo systems (only OS Galileo signals are used in this investigation, i.e., L1, E5a, and E5b) and a combined multiple-frequency GPS and Galileo system. It is tested with simulated data as real GPS L5 (only available in one Block IIR-M satellite now) and Galileo data are not yet available. The results, possibilities and weaknesses of the method are analysed. Note that the results of this investigation are not affected by the movement of the roving receiver because single-epoch data processing is used. Moreover, ambiguity resolution is not the main research interest in this paper, which assumes that ambiguities are fixed before the proposed algorithms are applied therefore no ambiguities are simulated in the carrier-phase data.

Cocktail multiple outlier detection algorithm

For high precision GPS positioning, the observation equations used in parameter estimation are usually linearized as:
E( y) = l + v = Ax (1)
E(v) = 0, l D (v) =D (l) =σ2 Cl (2)
where E( ) denotes the expectation operator, D( ) denotes the dispersion operator, l denotes the vector of double difference carrier phase observations, v denotes the vector of “true” residuals, A denotes the Jacobian matrix, x denotes the vector of unknown parameters, denotes the measurement variance factor, and denotes the cofactor matrix of observations. The unbiased estimate xˆ of x is based on double-difference single-epoch least squares solutions. Such a solution has known statistical properties and an expectation of the residuals of zero. It is given by the following formula.
xˆ = (AT WA)-1 ATWl
with
W = Cl -1     (3)
or in simplified form:
xˆ = N-1b with N = (ATWA),b = ATWl    (4)
where W, the weight matrix, is the inverse of the covariance matrix of observations. Observations contaminated by un-modelled error are generally referred to as outliers and will have the following statistical characteristics:
E=(vl) # 0 and/or D(vi )> σ2 Cl,ii    (5) Since a multipath error is an unmodelled error in the observation model and since more observations will be available in the modernised GPS and/or the new Galileo, the hypothesis that is to be tested here is that the increased redundancy might make multipath detection possible by statistical means. The author proposes to consider a phase observation affected by multipath as an outlier and seek to detect it by statistical testing. Cross (1994), based on Baarda (1968), which first introduced the classical outlier detection method widely used in geodetic network analysis, describes the test statistic ŵi for the ith correlated observation (i.e. where the weight matrix W contains off-diagonal elements) as:
(6)
where Cvˆ = W-1_ AN-1 AT    (7)
ei = [0 0 … 1 … 0 0]T    (8)
vˆ = Axˆ _ l    (9)
The test statistic is tested against a critical threshold c to test if contains a single multipath error whilst all other observations have only random and normally distributed errors
εi : H0 :li = lˉi +εi    (10)
H1: l = lˉ i +εi + Δi    (11)
where lˉ i denotes the true value of the quantities that have been observed. The Tau rejection criterion c described in (Pope, 1976) is used as the critical threshold and the observation is rejected when:

with
c = f (TN, TU,α )    (13) where NT is the number of (non-outlying) observations, NU is the degree of freedom, and α is the desired probability of Type I error. By shifting the position of the unity in the vector ei of Eq. 8, this method can be used to detect outliers in all measurements. The author named this method as MOD in this paper. Investigations described in (Lau, 2005) show that MOD is not sufficiently robust to tackle frequency-dependent multipath errors because it cannot handle the worst case scenario when multiple frequency multipath errors from a particular satellite are all (or any two of the three frequencies) in-phase. Therefore a Cocktail Multiple Outlier Detection (CMOD) algorithm is proposed in this section. The basic principle of this algorithm is based on the MOD method. In the MOD method, the test statistic of each measurement is tested against a critical threshold. However, the CMOD algorithm simultaneously tests all residuals of each satellite against a critical threshold and performs the test for all satellites in sequence. The test statistic (absolute value) of the measurements made on the three frequencies φf1 , φf2 and φf3 of a particular satellite s from a GNSS system is given by:
(14)
where Cvˆ = W-1_ AN-1 AT (7) es = [1 0 … 1 0 … 1 0 … ]T    (15)
in which the measurements of n satellites are arranged as:
[φs1 f1 … φsn f1 φs1 f2 … φsn f2 φs1 f3 … φsn f3 ]T    (16)
and vˆ = Axˆ _ l (9) Although Eq. 14 is theoretically restricted to the detection of only one outlier in the measurements, here the author tests the use of the formula to detect more than one. The author appreciates that this is theoretically incorrect, which is why the author uses the word ‘practical’ in the title of the paper. In order to tackle all possible cases when multipath errors of two or more frequencies from a satellite are in-phase as described in the previous section and shown in Fig. 1, the CMOD algorithm carries out statistical tests for residuals from all possible combinations of frequencies of a satellite. This is done by using different combinations of frequencies in the vector es(with the measurements are arranged as in Eq. 16) as: Test I: This test aims at the detection of the measurements of a satellite contaminated by significant in-phase multipath errors in all the three frequencies at the same time, as illustrated in the two black dotted lines of Fig. 1. es = [1 0 … 1 0 … 1 0 …]T (17) Test II: This test is for significant inphase multipath errors on both f1 and f2. es = [1 0 … 1 0 … 0 0 …]T (18) Test III: This test is for significant inphase multipath errors on both f1 and f3. es = [1 0 … 0 0 … 1 0 …]T (19) Test IV: This test is for significant inphase multipath errors on both f2 and f3. es= [0 0 … 1 0 … 1 0 …]T (20) Tests II to IV are intended to detect and reject serious multipath errors in any two of the three frequencies’ data of a satellite in cases where the multipath error in another frequency is insignificant. An example is shown in the orange dotted line of Fig. 1. The final test is required to detect any multipath error from just one frequency as indicated with the brown line in Fig. 1. This vector ei is the same as the one used in the MOD method: Test V: This test is for significant multipath errors on (f1, f2, or f3). es = [1 0 … 0 0 … 0 0 …]T (8) Tests I to V are performed for all combinations of all satellites at every epoch in order to detect multipath errors on all combinations of frequencies.

Description of testing

Description of simulated test datasets

The first two simulated datasets are referenced to the setup of a real multipath experiment carried out at the Laboratoire Central des Ponts et Chaussées (LCPC) near Nantes in France during May 2002 [6]. The real data was used in the validation of the multipath model used in the simulation of multipath data for the investigations of this paper, the validation results show that the model can generate very realistic phase multipath error (Lau and Cross, 2007). In the experiment, two Leica Geosystems System 530 receivers attached to lightweight AT502 antennas were used and a 5m by 2.5m steel panel

Fig. 2 Diagram showing the geometry of the LCPC test datasets

was constructed and placed about 5m to one of the receiving antennas in order to create a sufficiently large multipath signal. The antenna-reflector geometry is shown in Fig. 2. The geometry of the first two simulated datasets are the same but the reflectors with the assumed relative permittivities of 10 (e.g. flint glass) and 20 are used in multipath simulation. Therefore, the test datasets are denoted as LCPC-10 and LCPC-20 according to the relative permittivity used. The baseline length is about 86 m. The information used in multipath simulation is summarised in Table 1. The sky plot of the satellites-reflector-antenna geometry is shown in Fig. 3. Moreover, details of the satellites whose signals are contaminated by multipath are shown in Table 2.

Another virtual test site is at the IGS global tracking station LBCH in Long Beach, United States. In the LBCH test dataset, a concrete wall (reflector) is set 5m to the north of the antenna and a reference station 100m to its south. A relative permittivity of 7 (see Table 1) was used in the simulation of multipath in order to create strong multipath. This value is based on Stavrou and Saunders (2003), which found that the real part of the complex permittivity of concrete varied from 6.2 to 7 (for signals in the range of 1 to 95.9 GHz). Table 2 summarises the multipath simulation.

Also a kinematic data set, denoted as K-HK7-300, was simulated based on a 720m trajectory along a real street in Hong Kong. The roving antenna was assumed to set about 29cm (measured to L1 phase centre) above a 1m by 1m steel carry platform that travel at 1ms-1. The buildings were assumed to be made of concrete with a relative permittivity of 7. The reference station was assumed to be about 500 m to the north of the road. Table 3 summarises the simulation – note that some satellites are blocked for some of the time.

In the all of the test datasets, Leica System 530 receivers and AT502 antennas were assumed to have been used in data collection. Moreover, normal distributed random measurement noise with the standard deviation of 1 mm was generated in each

Fig. 3 Sky plot of Galileo (underlined) and GPS satellites in the LCPC datasets; the green area represents the refl ector; the grey area indicates that the data from the satellites in this area are contaminated by multipath

phase measurement. It should be pointed out that it was necessary to assume the use of a particular receiver and antenna because factors such as receiver correlator spacing and antenna gain pattern have an impact on the simulated multipath. The author does not believe the choices the author has made affect the overall conclusions of this research. Testing methodology Five scenarios with different GNSSs or combinations of frequencies have been tested and results compared with the known positions in all cases:

• Scenario 1: the current single

frequency GPS data;

• Scenario 2: the modernised

dual-frequency GPS data,

• Scenario 3: the future threefrequency

GPS data,

• Scenario 4: the future OS threefrequency

Galileo data, and

Table 1 Information for multipath simulation in the static test datasets

Table 2 Multipathing satellite information in the static test datasets

Table 3 Multipathing satellite information in the kinematic test datasets

• Scenario 5: the future OS Galileo
+ GPS multiple-frequency data.

Note that the reason for Scenario 1 is that
many surveying and geodetic applications only use L1 data in fixed solutions. Although L2 may have been used at the ambiguity fixing stage it may be too noisy to contribute usefully to the final solution.