Hostname: page-component-74656f9884-skzmx Total loading time: 0 Render date: 2026-11-06T04:16:51.667Z Has data issue: false hasContentIssue false
  • English
  • Français

Multipath Mitigation by Global Robust Estimation

Published online by Cambridge University Press:  26 June 2008

Sergio Baselga  and
Luis García-Asenjo
Sergio Baselga
Affiliation:
(Polytechnic University of Valencia, Spain)
Luis García-Asenjo
Affiliation:
(Polytechnic University of Valencia, Spain)
Rights & Permissions [Opens in a new window]

Abstract

In the last decades, robust estimation has been widely applied to overcome the presence of gross errors in observations. Recently it has been shown that robust estimation, if computed appropriately, is able to cope with a larger amount of gross errors and also capable of providing reliable estimations in situations where systematic errors are present. It is the purpose of this paper to propose the use of global robust estimation for computing baselines strongly affected by multipath effects.

Keywords

robust estimationleast squaresmultipathglobal optimization

Information

Type
Research Article
Information
The Journal of Navigation , Volume 61 , Issue 3 , July 2008 , pp. 385 - 392
Copyright
Copyright © The Royal Institute of Navigation 2008

1. INTRODUCTION

Signal multipath is the result of signal reflections from objects in the environment (Braash Reference Braash, Parkinson and Spilker1996). This effect degrades the signal and may be able to impede ambiguity resolution. Consequently, many efforts have been devoted to avoid multipath, the most obvious being the selection of a location free from potential reflections and antenna design improvements.

However, for many practical situations, especially in surveying engineering works, multipath effects remain after observation thus entering into the estimation stage. Many multipath mitigation techniques have been developed to act in the estimation process, to cite just a few: analysis of signal-to-noise ratio (Axelrad et al Reference Axelrad, Comp and Macdoran1996), wavelet analysis (Souza and Monico Reference Souza and Monico2004, Satirapod and Rizos Reference Satirapod and Rizos2005, Zhong et al Reference Zhong, Ding, Zheng, Chen and Huang2007, Aram et al Reference Aram, El-Rabbany, Krishnan and Anpalagan2007), modelling by means of harmonic functions (Amiri-Simkooei and Tiberius Reference Amiri-Simkooei and Tiberius2007), multipath reconstruction (Lau and Cross Reference Lau and Cross2007) and a very long etcetera.

On the other hand, robust estimation foundations were developed by Andrews et al (Reference Andrews, Bickel, Hampel, Huber, Rogers and Tukey1972), Huber (Reference Huber1981) and others, and then widely applied to geodesy and surveying (Fuchs Reference Fuchs1982, Harvey Reference Harvey1993, Yang et al Reference Yang, Wen, Xiong and Yang1999, and Qingming and Jinshan Reference Qingming and Jinshan1999, for instance) to provide a reliable estimation in an adjustment process in which gross errors participate. In particular, robust estimation was successfully applied to GPS positioning in a pioneer work by Wieser and Brunner (Reference Wieser and Brunner2002). Their conclusion was that “if the data contain outliers the result of robust estimation is far better than that of standard least squares. Otherwise the results are equal.” However, they computed the estimator by the usual means of an iteratively reweighted least squares (IRLS) algorithm. Recently Baselga (Reference Baselga2007) has shown that the process of computing a robust estimator by means of an IRLS algorithm highly undermines its estimation capabilities: IRLS robust estimation turns out to be only a local optimization around the initial least squares solution. Since this initial solution may be highly contaminated by the presence of undesirable errors and the correct solution may be far from the initial one, a process of global optimization was proposed to compute the robust estimator. Furthermore, this Global Robust Estimation (GRE) scheme proved to be capable of avoiding not only gross errors but also systematic errors, and it was successfully applied to the estimation of single frequency GPS baselines affected by a strong ionospheric delay which were unsolvable by classic methods (Baselga and García-Asenjo Reference Baselga and García-Asenjo2008).

This paper will show how baseline solutions free from multipath influences can be obtained by GRE. First we will present the adjustment problem, i.e. the mathematical model to be considered along with a suitable minimizing function (the robust estimator), then the computing strategy to be applied will be presented and finally we will consider an illustrative example in which GRE succeeds whereas classic methods fail in avoiding multipath effects.

2. MATHEMATICAL MODEL

Let us restrict the formulation, just for simplicity, to the case of single frequency phase observations (L1 phases). Note, however, that the use of additional observation types would not modify the essentials of the estimation process that we are proposing.

Let us consider observed single phases of receivers i and j from satellites r and s. The double differenced phase observing equation may be written as

(1)
\nabla \rmDelta \phi _{ij}^{rs} \equals \nabla \rmDelta \rho _{ij}^{rs} \minus \lambda \nabla \rmDelta N_{ij}^{rs} \plus \nabla \rmDelta T_{ij}^{rs} \minus \nabla \rmDelta I_{ij}^{rs} \plus \varepsilon _{ij}^{rs}

where subscripts denote receivers and superscripts satellites, \nabla \rmDelta \phi _{ij}^{rs} is the observed double differenced phase, \nabla \rmDelta \rho _{ij}^{rs} the double differenced geometric range, \nabla \rmDelta N_{ij}^{rs} the double differenced ambiguity for wavelength λ and \nabla \rmDelta T_{ij}^{rs} and \nabla \rmDelta I_{ij}^{rs} respectively the double differenced tropospheric delay, which is usually modelled, and the double differenced ionospheric delay, which is usually neglected for short baselines (or which cancels out in the dual frequency case). Residuals \varepsilon _{ij}^{rs} are assumed to follow a Gaussian distribution centred on zero, which happens if no significant systematic errors are present in the functional model, i.e. if Equation (1) accurately represents the observation reality.

Fixing Cartesian geocentric coordinates for one end of the baseline (say X i, Y i, Z i), using any approximate coordinates for the other end of the baseline (X j0, Y j0, Z j0) and reformulating the problem as the determination of the coordinate differentials for the end of the baseline (dX j, dY j, dZ j) along with the double differenced ambiguities (\nabla \rmDelta N _{ij}^{rs}), one can express the above equation, after adding the modelled tropospheric delay to \nabla \rmDelta \rho _{ij\setnum{0}}^{rs} and reordering, as

(2)
\openup2\eqalign{\tab\left( {\left[ {{{\partial \rho _{j}^{r} } \over {\partial X_{j} }}} \right]_{\setnum{0}} \minus \left[ {{{\partial \rho _{j}^{s} } \over {\partial X_{j} }}} \right]_{\setnum{0}} } \right)dX_{j} \plus \left( {\left[ {{{\partial \rho _{j}^{r} } \over {\partial Y_{j} }}} \right]_{\setnum{0}} \minus \left[ {{{\partial \rho _{j}^{s} } \over {\partial Y_{j} }}} \right]_{\setnum{0}} } \right)dY_{j} \plus \left( {\left[ {{{\partial \rho _{j}^{r} } \over {\partial Z_{j} }}} \right]_{\setnum{0}} \minus \left[ {{{\partial \rho _{j}^{s} } \over {\partial Z_{j} }}} \right]_{\setnum{0}} } \right)dZ_{j}\cr \tab \quad \minus \lambda \nabla \rmDelta N_{ij}^{rs} \equals \nabla \rmDelta \phi _{ij}^{rs} \minus \nabla \rmDelta \rho _{ij\setnum{0}}^{rs} \minus \varepsilon _{ij}^{rs}}

where subscript 0 denotes particularization for the approximate coordinates (X j0, Y j0, Z j0).

Therefore, the system of double differenced phase equations can be expressed, as usual, by

(3)
{\bf Ax} \equals {\bf l} \plus \nu

where A is the coefficient matrix, x the vector of unknowns – which contains differential corrections to approximate coordinates (float values) and ambiguities (integer values in theory) – l is the vector of observed minus approximate values and finally ν is the vector of residuals.

This system of equations has to be considered along with weight matrix P. A reasonable choice may be taken, for instance, from the SIGMA-ε weight model (Hartinger and Brunner, Reference Hartinger and Brunner1999). After adjustment, baseline components will be obtained by means of final coordinate differences: ΔX=X j0+dX j−X i, ΔY=Y j0+dY j−Y i, ΔZ=Z j0+dZ j−Z i.

3. ESTIMATION PROCESS

We are in presence of a dual estimation process, that is, some unknowns (coordinates) are float numbers whereas others (ambiguities) have to be integers. Classically the solution can be obtained by several different methods but all of them involve the application, under one or other scheme, of the least squares estimator, which is highly vulnerable to gross and systematic errors. As mentioned, Wieser and Brunner (Reference Wieser and Brunner2002) proposed the use of robust estimators computed, as usual, by an IRLS scheme. We propose an alternative that according to Baselga (Reference Baselga2007) means a significant improvement: the use of global robust estimation to overcome the possibility that robust estimation by IRLS falls only in a local optimum in the most complicated cases.

3.1. Least squares estimation

The classic least squares estimation process, see e.g. Teunissen Reference Teunissen1995, is performed in three steps:

  • • Least squares solution of the system, Equation (3). This yields real-valued ambiguities, which result in an approximate estimation for the baseline coordinates, namely the “float solution”.

  • • Integer ambiguity determination. The search for the best set of integer ambiguities is performed in the neighbourhood of the least squares solution found in the previous step.

  • • Constrained least squares adjustment of the system with the integer ambiguities found in the previous step. The result is known as the baseline “fixed solution”.

Note that if prior assumptions, i.e. the validity of the model and the presence of only random errors in the observations do not hold, then least squares adjustment no longer yields the most likely solution. Such is the case, for example, in the presence of multipath effects. Robust estimation is then a powerful alternative.

3.2. Robust estimation

Robust estimation is based on minimizing a function – named estimator – in order to attain a solution maximally resistant to the non-fulfilment of the mathematical model or the appearance of gross or systematic errors in observations, so that it also yields a solution almost identical to that of classic least squares when the model is correct and the observations are affected only by random errors.

Let us apply one of the most successful robust estimators that is widely used in literature (Fuchs Reference Fuchs1982, Harvey Reference Harvey1993, etc.): the L1-norm. Its solution is the one that minimizes the sum of absolute values of the residuals

(4)
{\min}\ \sum {\left\vert {\nu _{i} } \right\vert}

Typically, the routine to compute such a problem is an IRLS scheme, that is, the estimator minimization (4) is performed by successive least squares adjustments computing for each iteration adapted weights based on the previous adjustment residuals. The equivalent weight function for the L1-norm estimator to be used within least squares is

(5)
\omega _{i} \equals {1 \over {\left\vert {\nu _{i} } \right\vert}}

Therefore, minimizing the sum of squared residuals with weights (5) yields the L1-norm minimization.

(6)
{\min}\ \sum {\nu _{i} } ^{\setnum{2}} \omega _{i} \equals {\min}\ \sum {\nu _{i} } ^{\setnum{2}} {1 \over {\left\vert {\nu _{i} } \right\vert}} \equals {\min}\ \sum {\left\vert {\nu _{i} } \right\vert}

However, the problem with robust estimation as IRLS is a computational one, for if the initial least squares solution lies far away from the correct solution then the iterative process may attain only a local optimum. It seems quite obvious that if one has decided to abandon least squares estimation in favour of robust estimation, then solving robust estimation by means of an iteratively re-weighted least squares process may not be the most desirable solution. In Baselga (Reference Baselga2007) it was shown that, in fact, robust estimation by means of an iteratively re-weighted least squares process provides only a sort of “local robust estimation”, whereas robust estimation's best capabilities are only achieved if it is dealt with as a global optimization problem, thus forgetting all reference to least squares. Below, we present a possible scheme for robust estimation as a global optimization problem, GRE in short.

3.3. Global robust estimation

Having decided to focus on the L1-norm estimator, the estimation problem can be stated as

(7)
\left\{ {\matrix{ {{\bf Ax} \equals {\bf l} \plus \nu } \cr {\min}\ \rmSigma \left\vert {\nu _{i} } \right\vert}} \right.

For convenience we have transformed the general system of correlated observations (3) into a system of uncorrelated observations (expressed here with the same notation) by means of singular value decomposition of the original weight matrix P and subsequent transformation of (3) by pre-multiplication with the eigenvectors matrix, a procedure through which the vector of unknowns x remains invariant. The resulting system of equations has also been reduced to unit weight.

We will now conduct a global optimization process in the search for the best vector x in the sense of least L1-norm. Several successful global optimization techniques have been developed up to the moment, among them Genetic Algorithms (GA) and the Simulated Annealing (SA) method. We will only sketch the algorithm for the use of SA, a method that tries to emulate the process of crystalline networks' self-construction provided a sufficiently slow decrease in temperature occurs. In any case, the reader is referred to Baselga (Reference Baselga2007) and Berné and Baselga (Reference Berné and Baselga2004) for further details. The computational scheme may be as follows:

  1. 1. Start with an initial random solution xowithin the search domain D.

    In the present case the search domain may be a sufficiently wide area centred in the initial least-squares float solution with boundaries x i±kσxi (for a sufficiently large k, 30 or 50 for instance). Note the different requirements for the solution vector elements: some need to be real values (coordinate increments) and the rest integer values (double differenced ambiguities).

  2. 2. Obtain an increment to the previous solution and thus a new solution (provided it is in the search domain) {\bf x}_{i \plus \setnum{1}} \equals {\bf x}_{i} \plus \rmDelta {\bf x}_{i}.

    Displacement Δxi has to be generated at random, usually from a Gauss distribution whose standard deviation σi decreases with iterations. For example, one may start with an initial σi0=kσxi/2 so that the search domain can be easily inspected and apply a cooling factor of 0·999999n, where n is the iteration number. The affordable cooling speed depends on the problem complexity. Pay attention again to the different requirements for the solution vector elements: some are real values (coordinate increments) whilst other have to be integers (ambiguities).

  3. 3. Accept the new solution if it yields a better value for the score function – in our case the robust estimator (4) – otherwise: accept it with probability p (a small figure such as 0·01) or reject it.

  4. 4. Return to step 2 until typical displacement sizes σxi are below the required precision (for example 0·001 m for coordinates).

The solution obtained is guaranteed to be the global optimum in probabilistic terms, i.e. if the cooling process has been sufficiently slow. At any rate, a sufficient condition for the global optimum attainment can be easily checked: repeated executions of the algorithm must provide the same result.

4. EXAMPLE

Three very short baselines forming a triangle were observed by GNSS techniques. One of its vertex, point A, was selected to be in an environment potentially affected by strong multipath, see Figure 1.

Figure 1.

Point A in a potential multipath environment.

The results in Table 1 were obtained after an observation time of one hour and a register interval of 30 seconds by means of the Trimble Geomatics Office 1.62 software. Multipath in point A is clearly manifest: the first and third baselines, those including A, have large reference variances and large standard errors, along with low ratios that make us doubt the validity of these determinations. However, the closure computation gives sensible results, which reinforces the validity of the baseline determinations.

Table 1.

Baseline solutions for an observation time of 1 hour by TGO.

Let us now reduce observation times to 8 minutes. An analogous computation by Trimble Geomatics Office yields the results in Table 2. In this case, the short observation period has prevented the computation procedure from obtaining a correct solution: the closure computation points towards some incorrect determination and the comparison with the solution obtained for a one hour observing period locates the incorrectness in baseline A-B. In addition, it is significant how the reference variance and the standard error have increased for this baseline in reducing the observation time. However, note that the ratio for this baseline has increased to a value of 4·2, which, deceivingly, could have suggested that this incorrect determination was trustworthy. Moreover, a careful inspection of the residuals shows that marginal values (e.g. these exceeding 2·5σ) happen many times more than expected for the Gaussian behaviour assumption.

Table 2.

Baseline solutions for an observation time of 8 minutes by TGO and component differences with respect to the 1 hour observation solution.

Finally, let us apply Global Robust Estimation to this example, in a similar fashion as we did in Baselga and García-Asenjo (Reference Baselga and García-Asenjo2008) for coping with ionospheric delays. In this case we have followed the Simulated Annealing method as sketched in Section 3.3. We are trying to benefit from the fact that robust estimation is maximally resistant towards the non-fulfilment of the functional model and the appearance of gross errors, thus exploiting its ability to obtain correct results from largely contaminated samples. In our case, the highly contaminated sample – in which the proportion of reflected signals (i.e. with a multipath term) with respect to direct signals is high – is, in addition, of a short observing time (8 minutes). Results and comparison with the one-hour observation computations are shown in Table 3.

Table 3.

Baseline solutions for an observation time of 8 minutes by GRE (authors software) and component differences with respect to the 1 hour observation solution.

As we can see, contrarily to the classic baseline estimation process, GRE succeeds in obtaining correct baseline components despite the multipath affected observations and the short observation period.

5. CONCLUSION

Robust estimation computed as a global optimization process, or Global Robust Estimation, is presented as a powerful alternative to deal with multipath affected observations. After revising robust estimation basics, both for the classic iteratively re-weighted least squares scheme and for the successful global optimization scheme, we proposed to solve the functional model of double-differenced phases jointly in one step considering the dual nature of unknowns (float-integer) by a global optimization method upon minimization of a robust estimator functional. An example showed GRE advantages versus the classic least squares determination.

References

REFERENCES

Amiri-Simkooei, AR, Tiberius, CCJM (2007). Assessing receiver noise using GPS short baseline time series. GPS Solutions, 11(1): 21–35.Google Scholar
Andrews, DF, Bickel, PJ, Hampel, FR, Huber, PJ, Rogers, WH, Tukey, JW (1972). Robust estimates of location.: Surveys and advances. Princeton University Press, Princeton, NJ.Google Scholar
Aram, M, El-Rabbany, A, Krishnan, S, Anpalagan, A (2007). Single frequency multipath mitigation based on wavelet analysis. The Journal of Navigation, 60, 281–290.CrossRefGoogle Scholar
Axelrad, P, Comp, CJ, Macdoran, PF (1996). SNR-based multipath error correction for GPS differential phase. IEEE Transactions on Aerospace and Electronic Systems, 32: 650–660.CrossRefGoogle Scholar
Baselga, S (2007). A Global Optimization Solution of Robust Estimation. Journal of Surveying Engineering, 133(3): 123–128.Google Scholar
Baselga, S, García-Asenjo, L (2008). GNSS differential positioning by robust estimation. Journal of Surveying Engineering 134(1): 21–25.Google Scholar
Berné, JL, Baselga, S (2004). First-order design of geodetic networks using the simulated annealing method. Journal of Geodesy, 78, 47–54.CrossRefGoogle Scholar
Braash, MS (1996) Multipath effects. In: Parkinson, BW, Spilker, JJ (eds) Global positioning system: theory and applications. American Institute of Aeronautics and Astronautics, Cambridge, vol 1, 547–568.Google Scholar
Fuchs, H (1982). Contribution to the adjustment by minimizing the sum of absolute residuals. Manuscripta Geodaetica, 7: 151–207.Google Scholar
Hartinger, H, Brunner, FK (1999). Variances of GPS phase observations: The SIGMA-ε model. GPS Solutions, 2(4): 35–43.Google Scholar
Harvey, BR (1993). Survey network adjustments by the L1 method. Australian Journal of Geodesy, Photogrammetry and Surveying, 59, 39–52.Google Scholar
Huber, PJ (1981). Robust statistics. Wiley, New York.Google Scholar
Lau, L, Cross, P (2007). Development and testing of a new ray-tracing approach to GNSS carrier-phase multipath modelling. Journal of Geodesy, 81(11): 713–732.CrossRefGoogle Scholar
Teunissen, PJG (1995). The least-squares ambiguity decorrelation adjustment: a method for fast GPS integer ambiguity estimation. Journal of Geodesy, 70(1–2), 65–82.Google Scholar
Qingming, G, Jinshan, L (1999). Generalized shrunken-type robust estimation. Journal of Surveying Engineering., 125(4), 177–184.CrossRefGoogle Scholar
Satirapod, C, Rizos, C (2005). Multipath mitigation by wavelet analysis for GPS base station applications. Survey Review, 38: 2–10.Google Scholar
Souza, EM, Monico, JFG (2004). Wavelet Shrinkage: high frequency multipath reduction from GPS relative positioning. GPS Solutions, 8(3): 152–159.Google Scholar
Wieser, A, Brunner, FK (2002). Short static GPS sessions: robust estimation results. GPS Solutions, 5(3): 70–79.Google Scholar
Yang, Y, Wen, Y, Xiong, J, Yang, J (1999). Robust estimation for a dynamic model of the sea surface. Survey Review, 35: 2–10.CrossRefGoogle Scholar
Zhong, P, Ding, XL, Zheng, DW, Chen, W, Huang, DF (2007). Adaptive wavelet transform based on cross validation method and its application to GPS multipath mitigation. GPS Solutions, 10.1007/s10291-007-0071-yGoogle Scholar
Figure 0

Figure 1. Point A in a potential multipath environment.

Figure 1

Table 1. Baseline solutions for an observation time of 1 hour by TGO.

Figure 2

Table 2. Baseline solutions for an observation time of 8 minutes by TGO and component differences with respect to the 1 hour observation solution.

Figure 3

Table 3. Baseline solutions for an observation time of 8 minutes by GRE (authors software) and component differences with respect to the 1 hour observation solution.