Single Frequency Multipath Mitigation Based On Wavelet Analysis

Mohammad Aram; Ahmed El-Rabbany; Sri Krishnan; Alagan Anpalagan

doi:10.1017/S0373463307004146

Single Frequency Multipath Mitigation Based On Wavelet Analysis

Published online by Cambridge University Press: 20 April 2007

Sri Krishnan and

Mohammad Aram: Affiliation:
(Ryerson University, Toronto) (Email: maram@ryerson.ca)
Ahmed El-Rabbany: Affiliation:
(Ryerson University, Toronto) (Email: maram@ryerson.ca)
Sri Krishnan: Affiliation:
(Ryerson University, Toronto) (Email: maram@ryerson.ca)
Alagan Anpalagan: Affiliation:
(Ryerson University, Toronto) (Email: maram@ryerson.ca)

Article contents

Abstract
INTRODUCTION
WHY WAVELET ANALYSIS
PROPOSED ALGORITHM
RESULTS
CONCLUSION
References

Rights & Permissions

Abstract

Multipath is still one of the major error sources that degrades the accuracy of GPS positioning. The amount of multipath is highly dependent on the antenna's environment, which makes it difficult to isolate. Usually there is at least one in-view satellite which is more susceptible to multipath, particularly the one with the lowest elevation angle. To increase the positioning the best satellites must be selected (i.e. by least square or multipath mitigation) for computing a position. In this paper we propose an algorithm which picks up the best satellites (when there are more than four satellites in view) based on wavelet analysis for calculating a position. In this experiment, code and carrier measurements were collected in 15-minute segments by exploiting a single frequency (L1), stationary, navigation-grade receiver in a high-multipath environment. The magnitudes of these pseudoranges were often inflated by multipath error. We then post-processed the received data by applying wavelet filtering to the residuals (code minus carrier) to approximate the multipath values, and compute the receiver's position based on the selected satellites. Satellites were selected based on the residual values. To compare the results with the raw measurements, statistical elements were computed. The results showed significant improvement in variance of the estimated positions and, most importantly, a normalization of the data scatter-distribution was observed.

Keywords

Multipath Wavelet analysis GPS

Type: Research Article
Information: The Journal of Navigation , Volume 60 , Issue 2 , May 2007 , pp. 281 - 290

DOI: https://doi.org/10.1017/S0373463307004146 [Opens in a new window]
Copyright: Copyright © The Royal Institute of Navigation 2007

1. INTRODUCTION

Of the many sources of GPS error, multipath is among the most difficult to mitigate due to a lack of physical models. Multipath disturbance is largely influenced by the receiver's environment since satellite signals can arrive at a GPS receiver via multiple paths due to reflections from nearby trees, terrain, buildings, and vehicles. While the emergence of signal improvements (e.g. GPS modernization) will increase the likelihood of more accurate relative positioning, multipath will remain a relatively important issue, making positioning less accurate in urban environments. Because multipath is dominated by the environment immediately adjacent to the receiver, we must seek receiver-autonomous solutions. The use of wavelet analysis on GPS signals has proved to be an effective way of mitigating multipath. Satirapod and Rizos [Reference Satirapod and Rizos1] applied a wavelet decomposition technique to extract multipath from GPS observations using data collected by three dual frequency receivers. Zhang and Bartone [Reference Zhang and Bartone2] applied a wavelet filter on dual frequency measurements as well. In this paper, a wavelet analysis technique is applied to data taken from a single stationary receiver operating on L1 frequency only. In order for a GPS receiver to determine its position, it has to receive time signals from four different satellites. The pseudorange, P(t), of the user from the four satellites can be determined with the help of signal transit times between the four satellites and the receiver [Reference El-Rabbany3],

(1)

$P \lpar t \rpar \equals c \lsqb {t_{r} \lpar t \rpar \minus t^{\setnum{3}} \lpar {t \minus \tau } \rpar} \rsqb \plus \varepsilon _{p}$

where c is the speed of the light, t _r(t) is the reception time, t is the GPS time (true time), τ is the total signal travel time, t ³(t−τ) is the transmission time from a satellite, and ε_p is the error added to the measurements. Considering the atmospheric effect including ionospheric delay d _ion, and tropospheric delay d _trop along with multipath delay dm _r³ on the received signals, Equation 1 can be broken down into more detail [Reference El-Rabbany3]:

(2)

$P \lpar t \rpar \equals \rho \lpar {t\comma t \minus \tau } \rpar \plus c \hskip1\lsqb {d_{r} \lpar t \rpar \plus d^{\hskip1\setnum{3}} \lpar {t \minus \tau } \rpar} \rsqb \plus c \lsqb {d_{r} \lpar t \rpar \minus d^{\hskip1\setnum{3}} \lpar {t \minus \tau } \rpar} \rsqb \plus d_{trop} \plus d_{ion} \plus dm_{r}^{\setnum{3}} \plus \varepsilon _{p}$

Here ρ is the geometric distance between the satellite and the observing point, d ³ and d _r are the satellite and receiver equipment delay respectively. Carrier phase measurements, Φ(t), can also be used to compute the pseudoranges, which are more accurate than code measurements (P(t)).

Equation 3 shows a representation of the phase measurements [Reference El-Rabbany3]:

(3)

$\eqalign{ \Phi \lpar t \rpar \equals \tab \rho \lpar {t\comma t \minus \tau } \rpar \plus c\hskip1 \lsqb {\delta _{r} \lpar t \rpar \plus \delta ^{\setnum{3}} \lpar {t \minus \tau } \rpar} \rsqb \plus d_{trop} \plus d_{ion} \plus c \hskip1\lsqb {dt_{r} \lpar t \rpar \plus dt^{\hskip1\setnum{3}} \lpar {t \minus \tau } \rpar} \rsqb \cr \tab \plus \lambda \lsqb {\phi _{r} \lpar {t_{\setnum{0}} } \rpar \minus \phi ^{\setnum{3}} \lpar {t_{\setnum{0}} } \rpar} \rsqb \plus \lambda N \plus \delta m_{r}^{\setnum{3}} \plus \varepsilon _{\phi } \cr}$

The resultant time error caused by the error sources creates inaccuracies in the measurement. Among the errors multipath creates an incremental time delay, thereby exaggerating the respective pseudorange P(t) [Reference Zogg4]. To estimate the value of multipath from Equation (2) another variable called residual is introduced. Pseudorange can be calculated from both pseudo random code (PRN) and carrier phase information. The difference between these two pseudoranges is called residual. These residuals are processed by the proposed algorithm to be de-noised and ultimately used to mitigate multipath error.

2. WHY WAVELET ANALYSIS

In signal processing we often work with a time series of data. To process the data we look into the low and high frequency contents. When working with stationary data a Fourier analysis is sufficient, whereas working with time varying data the Fourier analysis is not a suitable tool because it is not able to show the time localization. Fourier basis functions are localized in frequency but not in time. Small frequency changes in the Fourier transform will produce changes everywhere in the time domain. Wavelets are local in both frequency (scale) via dilations and in time via translation. This localization makes wavelet transform a suitable tool for time varying signals. It is well known that the computational complexity of the fast Fourier transform is O(nlog ₂(n)), while for wavelet transform the computational complexity goes down to O(n)[Reference Sidney Burrus and Gopinath5]. As opposed to sinusoidal waves, which are considered to be big waves of infinite duration, wavelets are small waves, finite in duration. Where sinusoids are smooth and predictable, wavelets tend to be irregular and asymmetric. As Fourier analysis breaks up a signal into sine waves of various frequencies, wavelet analysis breaks up a signal into shifted and scaled versions of the wavelet. Generally speaking there are two types of wavelets.

1. Continuous wavelet transform (CWT)
2. Discrete wavelet transform (DWT)

The Fourier transform extracts from the signal details of the frequency content but loses all information on the location of a particular frequency within the signal. Time localization must then be achieved by first windowing the signal, and then by taking its Fourier transform. The problem with windowing is that the slice of the signal that is extracted is always the same length. Thus, the time slice (number of data points) used to resolve a high frequency component is the same as the number used to resolve a low frequency component [Reference Williams and Amaratunga6]. In contrast to windowed Fourier transform, the wavelet adapts the width of its time-slice according to the frequency components being extracted so that the resolution is high for high frequency components and low for low frequency components, because a wavelet uses short windows at high frequency and long window at low frequency.

2.1. Continuous wavelet transform

The CWT can be thought of as an inner product of the original signal with scaled shifted versions of the basis wavelet function ψ(t) [Reference El-Shimy, Osman, Nassar and Noureldin7]:

(4)

$\langle {x\comma \psi _{a\comma b} } \rangle \equals CWT \lpar {a\comma b} \rpar \equals \int \hskip-1{x \lpar t \rpar\psi \lpar t \rpar dt}$

(5)

${\rm where\colon }\quad \psi _{a\comma b} \equals {1 \over {\sqrt a }}\psi \left( {{{t \minus b} \over a}} \right)\hfil$

In Equation 5a represents the scale (dilation) and b is the time-shift (translation) parameter. Therefore, the wavelet transform of a continuous (analogue) signal x(t) is known as CWT which is defined as:

(6)

$CWT \lpar {a\comma b} \rpar \equals {1 \over {\sqrt a }}\int_{ \minus \infty }^{\infty } {x \lpar t \rpar \psi \left( {{{t \minus b} \over a}} \right) dt}$

The wavelet functions have time-widths adapted to their frequency such that at high frequency the time widths are narrow, while those at low frequency are much wider.

2.2. Discrete Wavelet Transform

Most time series are sampled as discrete values. The scaling and shifting variables are discretized so that wavelet coefficients can be described by two integers, m and n [Reference El-Shimy, Osman, Nassar and Noureldin7]. Thus, the DWT is given as:

(7)

$DWT \lpar {m\comma n} \rpar \equals {1 \over {\sqrt {a_{\setnum{0}}^{m} } }}\sum\limits_{k} {x \lsqb k \rsqb} \psi \lsqb {a_{\setnum{0}}^{ \minus m} n \minus k} \rsqb$

where x[k] is a signal or a digitized version of an analogue signal with sample index k, and ψ[n] is the wavelet. With different choices of m we obtain a geometric scaling: 1, 1/a ₀, 1/a ₀², … . It is found in practice that the most convenient value of a ₀ is 2 [Reference Ahmadi8]. This analysis method then consists of decomposing a signal into components at several frequency levels that are related by powers of two (a dyadic scale) [Reference El-Shimy, Osman, Nassar and Noureldin7].

In wavelet analysis we often speak of approximation and details. The approximations are the high-scale, low frequency components of the signal. The details are the low-scale, high frequency components [9]. Figure 1 shows a typical data analysis by wavelet.

Figure 1. A typical block diagram of wavelet analysis.

2.3. Thresholding

In wavelet decomposition when details are small, they might be omitted without substantially affecting the original signal. Therefore, de-noising refers to manipulation of wavelet coefficients for noise reduction in which coefficient values below a carefully selected threshold level are replaced by zero after which an inverse transform of modified coefficients is used to recover the de-noised signal. Mathematically, thresholding can also be described by a transformation of the wavelet coefficients in which transform matrix is a diagonal matrix with elements 0 or 1 [Reference Ahmadi8]. Zero elements force the corresponding coefficient below a given threshold to be set to zero while the others correspond to one, thus reducing the coefficients by the given threshold. The two approaches considered for de-noising are hard and soft thresholding. In hard thresholding only those wavelet coefficients with absolute values below or at the threshold level are affected. They are replaced by zero and others are left unchanged

$Hard\ Thresholding \left\{\hskip-2 {\matrix{ {W_{m} \equals \lpar W \hskip1\rpar} \cr {W_{m} \equals 0 }\hfill \cr} } \right.\ \matrix{ {if \ \vert W \vert \ges th}\hfill \cr {if \ \vert W \vert \lt th} \cr}$

where W and W_m are the wavelet coefficient before and after thresholding respectively. In soft thresholding coefficients above threshold level are also modified where they are reduced by the amount of threshold. Note that Donoho [Reference Donoho10] refers to soft thresholding as shrinkage since it can be proved that reduction in coefficient amplitudes by soft thresholding, also results in a reduction of the signal level thus a shrinkage [Reference Ahmadi8].

$Soft\ Thresholding \left\{\hskip-2 {\matrix{ {W_{m} \equals sign \lpar W \rpar \lpar {\vert W \vert \minus th} \rpar} \cr {W_{m} \equals 0}\hfill \cr} } \right.\ \ \matrix{ {if \ \vert W \vert \ges th} \hfill\cr {if \ \vert W \vert \lt th}\hfill \cr}$

3. PROPOSED ALGORITHM

The input data (residuals) are passed through db7 wavelet filters (level 5) to approximate multipath values. Then the best satellites are selected based on the least value of the multipath. Calculated positions are then outliered by 3σ. Figure 2 shows a block diagram of the proposed algorithm.

Figure 2. Block diagram of the proposed algorithm.

Seven data sets for one of the downtown Toronto locations (Yonge and Gerrard) were collected in 15-minute segments (with no cycle slips) at various times of the day over the period of one week in order to investigate the consistency of filters on varying multipath circumstances. The input data to be analyzed is a set of residuals (i.e. code minus carrier) given by Equation 8. When such a difference is performed all common errors such as satellite clock error, tropospheric error, and receiver clock error are eliminated. What remains is predominantly twice the ionospheric error, the pseudorange multipath error, and the carrier phase ambiguity. The carrier phase multipath is ignored since its values are very small compared to pseudorange multipath.

(8)

$P \lpar t \rpar \minus \Phi \lpar t \rpar \equals 2d_{ion} \plus dm_{r}^{\setnum{3}} \minus \lambda N$

The ambiguity term could be seen as a bias value which is removed by taking the mean value of the observations and subtracting the mean from the data. Therefore, the remaining values represent the ionosphere plus the multipath errors. To approximate the multipath error the wavelet analysis tool is applied to the remaining sum.

3.1. Stage 1: Multipath Approximation

A set of residuals denoted by Res_raw, collected from location L4, is shown in Figure 3. The wavelet procedures start with passing the signal (sequences of Res_raw) through a digital lowpass filter with impulse response g(n) and a highpass filter h(n). Filtering a signal corresponds to the mathematical operation of convolution of the signal with the impulse response of the filter. The convolution operation in discrete time is defined as follows:

(9)

${Re} s_{raw} \ast g \lpar n \rpar \equals \sum\limits_{k \equals \minus \infty }^{\infty } \hskip-3{{Re} s_{raw} \cdot g \lpar {n \minus k} \rpar}$

The lowpass filtering removes the high frequency information, but leaves the scale unchanged. Only the sub-sampling process changes the scale. Resolution, on the other hand, is related to the amount of information in the signal, and therefore, it is affected by the filtering operations. The lowpass filtering halves the resolution since the high frequencies are filtered, but leaves the scale unchanged. The signal is then sub-sampled by two since half of the number of samples are redundant [Reference Sidney Burrus and Gopinath5]. This doubles the scale. This procedure can be expressed mathematically as:

(10)

${{Re}} s_{low} \lpar k \rpar \equals \sum\limits_{n} {g \lpar k \rpar} \cdot {{Re}} s_{raw} \lpar {2n \minus k} \rpar$

(11)

${Re} s_{high} \lpar k \rpar \equals \sum\limits_{n} {h \lpar k \rpar} \cdot {Re} s_{raw} \lpar {2n \minus k} \rpar$

where Res _low(k) and Res _high(k) are the outputs of the lowpass and highpass filters respectively after sub-sampling by a factor of two. The above procedure, which is also known as sub-band coding, can be repeated for further decomposition. At every level, the filtering and sub-sampling will result in half the number of samples (hence half the time resolution) and half the frequency band spanned (hence double the frequency resolution). The filter used in the algorithm is from the Daubechis series (db7) [Reference Daubechies14] and decomposition level is considered to be five.

Figure 3. Time domain series of residuals (PRN 10).

We noted that successive approximations became progressively less noisy. Optimal de-noising needed a more subtle approach like thresholding. Thresholding involved discarding only the portion of the details or approximations that exceed a certain limit. The threshold limit in this experiment was set to the standard deviation (σ) of the wavelet coefficient achieved from filtering the residuals. Both soft and hard thresholding [Reference Donoho10] were applied to the wavelet coefficients, and coefficients were then reconstructed. In addition to de-noising the input data, the output is the approximated value of the pseudorange increment due to multipath since the ionosphere error is a slowly varying process while here we are dealing with a time varying error as shown in Figure 4. In Figure 4 the raw residuals form the trace with more fluctuation whilst the smoother trace represents the approximated values of the multipath for the corresponding satellite during the specified time.

Figure 4. Multipath approximation (PRN 10).

To validate the multipath approximation from wavelet filtering, a few data sets were collected by a dual frequency receiver. Multipath error was then computed and plotted using the TEQC software [11] as shown in Figure 5. Then from the same dual frequency data but using only the L1 frequency information, the multipath error was estimated by wavelet filtering as illustrated in Figure 6. By looking at the absolute values of both the TEQC and wavelet filtering, results are consistent with each other. Since some of the data were missing while using the L1 frequency data, there is a time shifting difference between Figures 5(b) and 6(b). Figures 7(a) and 7(b) show the frequency domain of the raw data (residuals before filtering), and the de-noised data respectively. It can be observed that the wavelet filtering has successfully de-noised the residuals. Figures 8(a) and 8(b) are two other samples of residuals with lower amount of noise compared with Figure 7(a).

Figure 5. Approximated multipath from TEQC using dual frequency data.

Figure 6. Approximated multipath from wavelet filtering, using single frequency data.

Figure 7. Residual spectrum before and after filtering.

Figure 8. Residual spectrum before and after filtering.

3.2. Stage 2: Finalizing the Computed Positions

To compute a position from the code measurements, we set a criterion based on the multipath (single frequency) value. In case there were more than four satellites in view, only the best four were selected to compute the position. We selected the ones with the lower multipath values at each epoch. Then outlier detection was applied to the computed positions. In this experiment, the outlier detection was based on the standard deviations (σ) in Easting and Northing components. These were computed and position estimates exceeding 3σ from the mean were rejected. The output from the second stage filter is the final position estimate as shown in Figure 9.

Figure 9. L4 positions: Raw (red) and filtered (black) data.

4. RESULTS

Figure 9 shows two data sets from the Yonge-Gerrard location. In each figure, the red dots represent the positions computed from the raw pseudoranges (data before filtering) while the black dots represent the positions after filtering. Figs. 9(a) and 9(b) show the positions for two different collection times. The horizontal and vertical axes represent the Easting and Northing respectively after normalizing to the scatter mean. Since satellite constellation and environmental conditions, such as traffic, change with time, the positions shown in plots (a) and (b) of Figure 9 scattered differently. A sample of the statistical measures are shown in Tables 1 and 2.

Table 1. Statistical results taken from Figure 9(a).

Table 2. Statistical results taken from Figure 9(b).

5. CONCLUSION

We have developed a two-stage wavelet filtering tool to mitigate GPS multipath. The first stage de-noises residuals to estimate multipath errors that overstate the magnitude of L1 pseudorange values. The second stage detects the outliers for a final position estimate. The results from our experiment show that excluding satellites with higher multipath could increase the position accuracy. The post filtered position scatters show smaller variance, eigenvalues, and kurtosis. On visual inspection, our approach to wavelet-based filtering provides a strong degree of statistically measurable noise removal in the seven, 15-minute datasets studied. De-noising of signals from a stationary receiver has applications in low-cost survey and location applications.

ACKNOWLEDGEMENTS

The authors would like to thank Mr Bern Goush and Mr Aiden Morrison for their help during the initial phase of the project.

References

REFERENCES

[1]Satirapod, C., Rizos, C., (2003). “Multipath Mitigation by Wavelet Analysis for GPS Base Station Applications”, Survey Review, Vol. 38, No. 295.Google Scholar

[2]Zhang, Y., Bartone, C., (2004) “Real-time Multipath Mitigation with WaveSmooth Technique using Wavelets”, ION GNSS 17th International Technical Meeting of the Satellite Division, pg 1181–1193, 21–24 September 2004.Google Scholar

[3]El-Rabbany, A., (2006) “Advanced Satellite Positioning”, Unpublished Lecture Notes, Department of Civil Engineering, Ryerson University.Google Scholar

[4]Zogg, J., (2002) “GPS Basics, Introduction to the System Application Overview”, u-blox Corporation, www.telecom.fh-htwchur.ch/zogg/Dateien.Google Scholar

[5]Sidney Burrus, C., Gopinath, A. Ramesh, (1998) “Wavelet and Wavelet Transforms”, Rice University, Prentice-Hall, Inc.Google Scholar

[6]Williams, J. R. and Amaratunga, Kevin, (2001) “Introduction to Wavelets in Engineering”, Intelligent Engineering Systems Laboratory, Massachusetts Institute of Technology, Cambridge.Google Scholar

[7]El-Shimy, N., Osman, A., Nassar, S., Noureldin, A., (2003). “Wavelet Multiresolution Analysis”, GPS World.Google Scholar

[8]Ahmadi, H., (2004) “Wavelets: Principles and Applications in Signal Processing”, University of Tehran, Electrical and Computer Engineering Department, Lecture Notes.Google Scholar

[9]The MathWorks, Inc., http://www.mathworks.com/access/helpdesk/help/techdoc/matlab.shtml.Google Scholar

[10]Donoho, D. L., (1995) “De-noising by Soft-Thresholding”, IEEE Trans. on Information Theory, Volume 41, Issue 3.CrossRef Google Scholar

[11]UNAVCO, TEQC software, http://facility.unavco.org/software/teqc/teqc.html.Google Scholar

[12]Misra, P., Burke, B., Pratt, M., (1999) “GPS Performance in Navigation”, Proceedings of the IEEE, VOL. 87.Google Scholar

[13]Bisnath, S., Langley, R., (2001). “Pseudorange Multipath by Means of Multipath Monitoring and De-Weighting”, presented at KIS 2001, 5–8 June 2001, Alberta.Google Scholar

[14]Daubechies, I., (1992) “Ten Lectures on Wavelets”, CBMS-NSF Regional Conference Series in Applied Mathematics, Vol. 91, SIAM, Philade.Google Scholar