Wavelet packets based denoising method for measurement domain repeat-time multipath filtering in GPS static high-precision positioning

Repeatable satellite orbits can be used for multipath mitigation in GPS-based deformation monitoring and other high-precision GPS applications that involve continuous observation with static antennas. Multipath signals at a static station repeat when the GPS constellation repeats given the same site environment. Repeat-time multipath filtering techniques need noise reduction methods to remove the white noise in carrier phase measurement residuals in order to retrieve the carrier phase multipath corrections for the next day. We propose a generic and robust three-level wavelet packets based denoising method for repeat-time-based carrier phase multipath filtering in relative positioning; the method does not need tuning to work with different data sets. The proposed denoising method is tested rigorously and compared with two other denoising methods. Three rooftop data sets collected at the University of Nottingham Ningbo China and two data sets collected at three Southern California Integrated GPS Network high-rate stations are used in the performance assessment. Test results of the wavelet packets denoising method are compared with the results of the resistor–capacitor (RC) low-pass filter and the single-level discrete wavelet transform (DWT) denoising method. Multipath mitigation efficiency in carrier phase measurement domain is shown by spectrum analysis of two selected satellites in two data sets. The positioning performance of the repeat-time-based multipath filtering techniques is assessed. The results show that the performance of the three noise reduction techniques is about 1–46 % improvement on positioning accuracy when compared with no multipath filtering. The statistical results show that the wavelet packets based denoising method is always better than the RC filter by 2–4 %, and better than the DWT method by 6–15 %. These results suggest that the proposed wavelet packets based denoising method is better than both the DWT method and the relatively simple RC low-pass filter for noise reduction in multipath filtering. However, the wavelet packets based denoising method is not significantly better than the RC filter.


Introduction
Multipath is one of the most important Global Positioning System (GPS) error sources in high-precision positioning.Multipath errors are caused when direct signals from satellites are mixed with those reflected from objects in the vicinity of the antenna.Theoretically, the maximum amplitude of multipath error in a phase measurement is a quarter of the observing wavelength, for example, it is about 5 cm for the GPS L1 carrier.Therefore, multipath mitigation is crucial to achieve centimeter and millimeter positioning accuracy.Details of carrier phase multipath effect can be found in Lau andCross (2006b, 2007).
Multipath mitigation techniques can be classified into site-dependent, hardware-dependent, and algorithm-dependent techniques.Park et al. (2002) and Wanninger and May (2000) describe in situ multipath calibrations for reference stations based on the repeatable satellite-reflector-antenna geometry in about one sidereal day.Another site-dependent technique is multipath environmental modeling; this technique estimates multipath errors based on the known satellite-reflector-antenna geometry at static GPS stations.Lau and Cross (2007) describe a rigorous ray tracing algorithm using the known satellite-reflector-antenna geometry and physical properties of reflectors to determine phase multipath errors; this algorithm can be used for multipath mitigation at static GPS stations.Site-dependent techniques can mitigate multipath effectively since multipath errors can be determined from the known and/or repeatable satellite-reflector-antenna geometry.
Hardware-dependent multipath mitigation techniques can be categorized into antenna based and receiver based techniques.Antenna based techniques can be special antenna designs such as choke ring antennas (Filippov et al. 1998) and Trimble's Zephyr antennas (Krantz et al. 2001).Antenna gain patterns can also be used to mitigate multipath by reducing the gain of a low-elevation signal.This method is based on the increase of multipath as the satellite elevation angle decreases.However, such a relationship is not always true, especially for the phase multipath error for which the magnitude and phase change with varying differential path delay, i.e., the path difference between the direct and indirect signal paths changes with the satellitereflector-antenna geometry (Lau and Cross 2007).Ray et al. (2001) and Lau and Cross (2006a) describe antenna array techniques based on the geometric correlation of multipath errors at closely spaced antennas.Be ´taille et al. (2006) describe a receiver based technique, called the phase multipath mitigation window that relies on the gated correlator.This technique can effectively mitigate multipath with the differential path delay of more than 7.5 m.A similar technique, called vision correlator, is described in Fenton and Jones (2005).This correlator can effectively mitigate multipath with the differential path delay of more than 5 m.However, those two receiver based techniques may increase the noise level in measurements and reduce the signal-to-noise ratio (SNR).
Regarding algorithm-dependent techniques, Wieser and Brunner (2002) describe the SIGMA-F model which controls the parameter estimation by iteratively re-weighted least squares.The core of this method is the fuzzy system that uses the conventional outlier detection method and the measured SNR as the data quality indicator to control the re-weighting scheme.Lau and Cross (2006b) describe the modified SNR based stochastic model.This method applies the modified SNR to form the stochastic model.This model was developed based on the fact that phase multipath error and SNR are orthogonal (Lau and Cross 2006b) and the use of measured SNR may provide wrong information for the quality of phase measurement.The modified SNR is determined by using calibrated and empirical data.It can fix the orthogonality problem, which means that the modified SNR can reflect the data quality directly and appropriately.Comp and Axelrad (1996) and Lau and Cross (2005) investigate the use of SNR to estimate phase multipath errors by adaptive filters with spectrum analyses; however, this technique is also affected by the orthogonality between phase multipath error and SNR.
Owing to the unique condition of repeatable satellite geometry in about one sidereal day at continuous static antennas, repeat-time-based filter may be the best technique to mitigate multipath effect that remains in measurements even though the multipath mitigating antennas and receivers are used in data collection.Genrich and Bock (1992) use sidereal filtering to correct the multipath contaminated position errors in continuous monitoring of crustal deformation.The fixed sidereal day of 23 h and 56 min is used to correlate the positioning errors in consecutive days.Seeber et al. (1998) investigate the exact period of the sidereal repeatability for individual satellites using cross-correlation of double differences in time domain, correlation of elevation and azimuth time series, and computation from individual ephemeris.Different repeat times of 240-256 s in sidereal advancement for satellites are found in the literature, and results show that the repeat times of satellites determined by the three methods are in good agreement.Axelrad et al. (2005) evaluate three methods for estimation of the repeat time.The methods are as follows: the computation of orbital period from the semimajor axis and the correction to the mean motion given in the broadcast ephemeris; computation of the repeat time by interpolating post-processed GPS orbit solutions to the equator crossing on subsequent days; and determination of the actual repeat geometry for a selected location and then identifying the associated time shift.Agnew and Larson (2007) also investigate two methods of finding the repeat times of the GPS constellation.The first method is the use of the semimajor axis and the correction to the mean motion given in the broadcast ephemeris to compute the orbital periods, and the other is the aspect repeat time of the topocentric positions of the satellites.The investigations show that the determined repeat time for satellites by various methods agree usually within 3 s.They suggest that the best corrections for multipath will be obtained by processing data from each satellite separately, rather than working with a time series of positions.Lau (2012) shows that the position and measurement domain multipath filtering techniques perform similarly in practice because multipath errors do not change too much in the time difference between the actual sidereal day (in the measurement domain) and the mean sidereal day (in the position domain) computed from the sidereal days of all observing satellites with the orbit parameters in the broadcast ephemeris, and because the advanced receiver correlators can filter out long-delay multipath signals that lead to fast change in differential path delay.
Wavelet-based techniques have been widely used in signal processing, data compression, data analyses in engineering, sciences, and finance, etc. Huang et al. (2003) apply wavelets in dynamic deformation monitoring for high-rise buildings.Zhang and Bartone (2004) apply wavelet decomposition for real-time multipath mitigation in code minus carrier observable.Wavelet-based techniques have also been used in repeat-time-based multipath mitigation for static receivers.Satirapod and Rizos (2005) use Symlets wavelet transform-based multiresolution analysis to retrieve multipath signals and the first three levels of decompositions for multipath disturbance are investigated.Souza and Monico (2004) use wavelet shrinkage technique to retrieve multipath signals and six mother wavelets are compared.Since there are many mother wavelets such as Haar, Coiflet, Daubechies, and Symlet wavelets and some of them have various numbers of coefficients such as Coiflet's C6, C12, C18, C24, etc., where the numbers after C represent the number of coefficients, we do not attempt to compare the performance of mother wavelets on retrieving multipath errors from noisy measurement residuals in some selected data sets.The actual performance of multipath error retrieval for a given data set depends on various parameters of signal and receiver combination, including the signal-type modulation, code chipping rate, the pre-correlation bandwidth and filter characteristics, the number of received multipath signals, the relative power of multipath signals, the differential path delay (long-or short-delay multipath), chip spacing between correlators, and the type of discriminator and algorithm used for code and carrier tracking (Lau and Cross 2007).Therefore, validations with few data sets may not be sufficient to obtain the best mother wavelet for the repeat-time-based multipath mitigation technique in all multipath situations.This work aims to propose a generic wavelet packet decomposition technique to denoise the carrier phase residuals obtained in the first day (Day 1) in order to retrieve the multipath corrections for the carrier phase measurements in the next day (Day 2).The performance of the wavelet packets based denoising technique is compared with that of the resistor-capacitor (RC) low-pass filter and the single-level discrete 1-D wavelet transform (DWT)-based denoising method.The RC filter is based on the exponentially weighted moving average (Roberts 2004), and the single-level DWT and inverse discrete 1-D wavelet transform (IDWT) functions of MathWorks MATLAB TM (MathWorks 2015) are used in the DWT method.The semimajor axis and the correction to the mean motion given in the broadcast ephemeris are used to determine the repeat time of the satellite.This repeat-time determination requires very little computing load and is quick, and produces similar repeat time as the other methods described above (Axelrad et al. 2005;Seeber et al. 1998).Moreover, this repeat-time determination method has the least latency when comparing with those methods that require precise orbits and those approaches involving correlation of residual time series between days.Therefore, this multipath filtering algorithm can be applied for realtime deformation monitoring.The three noise reduction techniques for repeat-time-based multipath filter are tested with five data sets.Spectrum analyses of two selected satellites in data sets are used to show the filtering efficiency of wavelet packets denoising.The root-mean-square (RMS) positioning errors (northing, easting, height) of the multipath filtering with the three noise reduction methods are compared with those of the standard least-squares single-epoch solution.Improvement of the three noise reduction methods for multipath filtering on positioning accuracy is used to assess their performance.
In the literature, most wavelet-based techniques in GNSS are applied to position domain sidereal filtering and position domain time series applications.Most contributions propose wavelet-based techniques, but no comparison with other similar techniques can be found.This work may be the first study that proposes a wavelet packets denoising technique in measurement domain repeat-time-based multipath filtering and carries out a rigorous test with comparison with two advanced denoising methods.

Repeat-time-based Multipath Filtering
Repeat-time multipath filtering is also known as sidereal filtering.As mentioned above, sidereal filtering can be carried out in the position domain or measurement domain.The theory of the two sidereal filtering approaches for short/medium baselines and the basic formulas involved are described below.
Position domain sidereal filtering uses position residuals obtained in first day (Day 1) as the correction to the position solutions in the next day (Day 2) according to the repeat time.Putting this approach in formulas, the positioning solution of each epoch in Day 1 can be written as: where N; Ê; Ĥ denote the best estimated positioning components in northing, easting, and height, respectively, of the static station, N; E; H represent the known northing, easting, and height, respectively, of the static station, M N ; M E ; M H denote the multipath contaminated northing, easting, and height, respectively, and e N ; e E ; e H are the random errors in northing, easting, and height, respectively.Putting the known position vector of (1) to the lefthand side, it becomes: Positioning component residuals on the right-hand side contain multipath errors and random errors.Multipath contaminated northing, easting, and height components in time series can be obtained by time/frequency response analysis, that is, obtaining the multipath errors from the position residuals.Denoising the position residuals can also be used to obtain the multipath errors, that is, removing the random noise from the position residual.In the position domain sidereal filter, position residuals in Day 1 are the input signal for the time/frequency response analysis and denoising.The corrected position of each epoch in Day 2 is obtained by: Measurement domain multipath filtering uses the individual repeat time to apply the carrier phase measurement residuals obtained from the first day (Day 1) as the correction to the measurements in the next day (Day 2) for each satellite.
The dual-frequency phase observables given in Strang and Borre (1997): where U kl 1;ij ; U kl 2;ij , and U kl 5;ij are the double-difference phase observations between satellites k and l, and stations i and j for L1, L2, and L5 carriers, respectively, q kl ij denotes the double-difference geometric range, I kl ij denotes the doubledifference ionospheric effect, T kl ij denotes the double-difference tropospheric delay, N kl ij denotes the double-difference integer ambiguity, and e kl ij denotes the doubledifference measurement noise.
For short or medium baselines, the double-difference residuals consist of mainly multipath errors and small measurement noise.Since carrier phase multipath error is always less than a quarter of the observing carrier wavelength, the measurement residuals obtained in Day 1 do not need to consider the integer part, which is the carrier phase ambiguity

Denoising with wavelet packet decomposition
In most applications, one discrete wavelet transform (DWT) is not enough, and often it is necessary to do a complete wavelet packet decomposition (Jensen and la Cour-Harbo 2001).This means applying the DWT several times to various signals.Compared with the conventional short-time FFT (SFFT), wavelet packets have the advantages that they are well localized in both time and frequency, and the sidelobe energy leakage is much smaller than that of the SFFT (Zhang and Dill 1999).Moreover, an advantage of wavelet packets over the standard wavelet transform is observed in Paiva and Galva ˜o (2008) in that the wavelet packet technique performs an adaptive partitioning of the frequency axis.Wavelet packets represent a generalization of the method of multisolution decomposition and comprise the entire family of sub-band-coded (tree) decompositions.The inverse transform of this decomposition has a perfect reconstruction property.This property is essential for the proposed generic waveletbased denoising algorithm in retrieving the multipath signals from measurement residuals.
The proposed wavelet packet denoising technique consists of the following steps: 1. Wavelet transform: transform the double-difference phase measurement residuals (input signal) to the wavelet domain using Daubechies 4 (D4) transform, 2. Three-level wavelet packet decomposition: repeat Step 1 three times to decompose all the wavelet coefficients (high-pass and low-pass coefficients), see Fig. 3, 3. Denoising by thresholding: apply a threshold to remove the noise, and 4. Inverse wavelet transform: inverse-transform the denoised wavelet coefficients to the signal domain (multipath errors).The output is the ''noise-free'' double-difference phase multipath errors.
Details of the above steps are described as follows: 1. Wavelet transform: Daubechies orthogonal wavelets D2-D20 (even index numbers only) are commonly used.The index number refers to the number N of coefficients.Each wavelet has a number of zero moments or vanishing moments equal to half the number of coefficients (Mallat 2009).Wavelet packets generalize the compactly supported wavelets of Daubechies (Coifman et al. 1990;Wickerhauser 1994); therefore, the transform with two vanishing moments from the Daubechies family (i.e., D4) is adopted in the proposed algorithm.Moreover, since the full wavelet packet decomposition and reconstruction are used in Steps 2 and 4, the selection of a mother wavelet becomes less important (Wickerhauser 1994).2. Three-level wavelet packet decomposition: Wavelet packets are used to gain the advantage of better frequency resolution representation.Three-level decomposition shown in Fig. 3 is selected in the proposed algorithm with a good balance between frequency resolution and processing load.The advantage of this further series of decompositions is that the time frequency plane is partitioned more precisely (Samantaray and Dash 2007).Since full wave packets use all the wavelet coefficients, it has a perfect reconstruction property for the inverse wavelet transform.This property is essential for the generic wavelet-based multipath filter proposed in this work.
In this step, we use 512 samples and partition complete where otherwise: The hard thresholding method is selected for the proposed algorithm to avoid changing the magnitude of ''multipath signal'' in the wavelet coefficients and because Souza and Monico (2004) show that hard thresholding is more suitable for GPS applications.Donoho and Johnstone (1994) have investigated thresholding in statistical approaches, and the universal threshold is proposed as: where n denotes the number of samples in the Gaussian distributed noise vector N and The noise vector in Level 1 decomposition is used for the threshold determination in the proposed algorithm.
All the Level 3 wavelet coefficients are denoised by the determined threshold.

Test data description
In order to assess the impact of repeat-time-based filtering on multipath mitigation clearly, short baselines are used in the tests.Three short baselines are used.The first one is a very short 10.5 m baseline on the roof of the Science and Engineering Building at the University of Nottingham Ningbo China (UNNC); the GNSS reference stations and the site environment are shown in Fig. 6.A Leica AR20 3D GNSS choke ring antenna is set on the north pillar, and it is connected to a Leica GR25 GNSS receiver.On the south pillar, a Leica AR25 GNSS choke ring antenna is connected to a Leica GR10 GNSS receiver.The north station acts as the reference station in the data sets of the UNNC baseline, a 10°satellite elevation mask, and 1 Hz data rate were used in the data processing.Three sessions of data were selected in this test, and each data set (session) was observed in two consecutive days.Details of the data sets are listed in Table 1.Dual-frequency double-difference phase observable was used in the data processing, and the highest elevation satellite was selected as the reference satellite for double difference.
Another two baselines are formed between three Southern California Integrated GPS Network (SCIGN)/ California Real Time Network (CRTN) high-rate (1 Hz data) stations in Los Angeles.The stations are LAND, MIDA, and POMM.Trimble NetRS GPS receivers and Ashtech Dorne Margolin with choke rings (ASH701945B_M SCIS) antennas are used at the three stations.POMM was selected as the reference station in the data processing, and the POMM-MIDA baseline length is about 1.8 km and about 2.6 km for the POMM-LAND baseline.A 10°satellite elevation mask and 1 Hz data rate were used in the data processing.The observation times of these data sets are from 01:00:00 to 03:00:00 (UTC) of January 1, 2015, and from 00:55:00 to 02:56:00 (UTC) of January 2, 2015.
The standard single-epoch least squares with predetermined carrier phase ambiguities is used in data processing.The repeat-time-based multipath filtering with the three noise reduction methods applies to all satellites in the data sets (i.e., no identification of multipathing satellite is needed).Performance assessment of the three noise reduction

Results and analyses
Time fast Fourier transform (TFFT) analysis on the original, i.e., no multipath filtering applied, and the wavelet packets denoising based multipath-filtered double-difference phase measurements of the satellite PRN02 in the UNNC Data set 1 and the satellite PRN09 in the POMM-LAND data set has been carried out.Moreover, spectrograms using short-time Fourier transform, a segment length of 256 samples in a Hamming window and an overlap of 167 samples, with MATLAB TM on the original (no multipath filtering) and the wavelet packets denoising based multipath-filtered double-difference phase measurements of the satellite PRN02 in the UNNC Data set 1 are shown in Fig. 9.All the above spectrum analysis results show that no significant multipath errors remain in the measurement residuals after the multipath filtering.
The RMS positioning errors of the UNNC Data sets 1, 2, and 3 are shown in Tables 2, 3, and 4, respectively.Percentage changes of the RMS errors of the multipath filtering with RC filter, the DWT method, and the wavelet packets based denoising are compared with the original least-squares solutions (no multipath corrections) are also shown in the tables.Moreover, the results for the POMM-MIDA and POMM-LAND baselines are shown in Tables 5  and 6, respectively.Regarding the computational load in data processing, both the RC filter and wavelet packets based denoising technique, coded in C?? by the author, are very efficient.The DWT method is coded in MATLAB TM , and it is very efficient as well.With the current computer technology, there is no noticeable difference in processing time between the three noise reduction methods in practice.
From the above tables, the three filters show 3D positioning accuracy improvements when compared with the  The results may be controversial in that the DWTbased sidereal filter does show positioning accuracy improvement but it is not better than the RC and wavelet packets denoising based sidereal filters.The overall performance of the wavelet packets denoising is always better than those of the RC-and DWT-based filters, but the wavelet packets filter is only slightly better than the RC filter (\5 %).However, only the wavelet packets based filter shows improvements in all positioning components, i.e., northing, easting, and height, in all the data sets.Gokhale and Khanduja (2010) assess the performance of the discrete wavelet transform (DWT) and the wavelet packet decomposition (WPD) and conclude ''The performance of wavelet packet is appreciable while comparing with the discrete wavelet transform decomposition technique since wavelet packet analysis can provide a more precise frequency resolution than the wavelet analysis.It also has compact support in time as well as in frequency domain and adapts its support locally to the signal which is important in time varying signal''.GPS multipath errors vary with time because of the change in satellite-reflector-antenna geometry.The finding in this work agrees with the finding in Gokhale and Khanduja (2010) on the DWT and wavelet packets performance.and wavelet packets denoising results are 16.6, 23.9, and 26.8 % respectively when compared with the original no multipath filtering results; the improvement is subject to the severity of multipath errors in the data sets/sites.The test results demonstrate that all filters significantly mitigate multipath errors and improve positioning accuracy, with the wavelet packets based multipath filter being always better than the other two methods.The results of this investigation suggest that the proposed wavelet packets based denoising method is better than relatively simple low-pass filters for noise reduction in multipath filtering although it is not significantly better (\5 % in all cases).Moreover, the results show that the wavelet packets based method is better than the DWT-based method in the repeattime-based multipath filtering.This is the case because the wavelet packets method performs better for time varying input signals, i.e., when amplitude and frequency of the input signal change with time; the amplitude and frequency of multipath errors change with time as the differential path delay and the satellite-reflector-antenna geometry change.Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/),which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
where _ N; _ E; _ H denote the multipath-corrected positioning components in northing, easting, and height.Since this technique applies corrections in the position domain, the different repeat time for satellites cannot be used.The mean repeat time, therefore, must be used in this technique, and it is obtained by taking the mean of all computed repeat times of all satellites in view at each observation epoch.The basic steps of the position domain sidereal filtering technique are shown in Fig.1.Examples of applying wavelets in position domain sidereal filters can be found inYe et al. (2013) andKhelifa et al. (2011).
. The basic steps of the measurement domain sidereal filtering technique are shown in Fig. 2. In theory, repeat-time-based multipath filtering increases the noise level in the position solutions of Day 2. Low-pass filter or other techniques must be used to filter out the noise in position and measurement residuals in the two domains of repeat-time-based multipath filter.Measurement domain sidereal filtering is the focus of this work.Comparison between the position domain and measurement domain sidereal filtering techniques can be found in Lau (2012).

Fig. 1
Fig. 1 Flowchart for the position domain sidereal filtering technique

Fig. 2
Fig. 2 Flowchart for the measurement domain sidereal filtering technique

Fig. 3
Fig. 3 Three-level wavelet packet decomposition used in the denoising algorithm.The input signal is decomposed three times to obtain all highpass and low-pass coefficients

Fig. 5
Fig. 5 Flowchart for the proposed wavelet packet denoising algorithm

Fig. 6
Fig. 6 GNSS reference stations in UNNC and the site environment

Fig. 7
Fig. 7 Positioning errors in northing (top), easting (middle), and height (bottom) show the multipath signal sinusoidal signature and noise of the UNNC Data set 1

Fig. 8
Fig. 8 Positioning errors in northing (top), easting (middle), and height (bottom) show the multipath signal sinusoidal signature and noise of the POMM-MIDA data set

Fig. 10
Fig. 10 TFFT analysis of Day 2 height errors in the UNNC Data set 1 (top), set 2 (middle), and set 3 (bottom).The unit in magnitude is millimeters

Fig. 11
Fig. 11 TFFT analysis of Day 2 height errors in the POMM-MIDA (top) and the POMM-LAND (bottom) data sets.The unit in magnitude is in millimeters

Table 1
Observation time of the test data sets for UNNC multipath data

Table 2
RMS positioning errors of the original least-squares solution, the RC filter-based multipath filter, the DWT method, and the wavelet packets based multipath filter for the UNNC Data set 1 Table3RMS positioning errors of the original least-squares solution, the RC filter-based multipath filter, the DWT method, and the wavelet packets based multipath filter for the UNNC Data set 2 Table4RMS positioning errors of the original least-squares solution, the RC filter-based multipath filter, the DWT method, and the wavelet packets based multipath filter for the UNNC Data set 3 Table5RMS positioning errors of the original least-squares solution, the RC filter-based multipath filter, the DWT method, and the wavelet packets based multipath filter for the POMM-MIDA baseline Table6RMS positioning errors of the original least-squares solution, the RC filter-based multipath filter, the DWT method, and wavelet packets based multipath filter for the POMM-LAND baseline