A New Method for Multipath Filtering in GPS Static High-Precision Positioning

It is well known that multipath is one of the main sources of errors in GPS static high precision positioning of short baselines. Most algorithms for reducing multipath manipulate the GPS double difference (DD) observation residuals as input signal in GPS signal processing. In the traditional multipath mitigation methods, applying the wavelet transform (WT) to decompose the GPS DD observation residuals for identifying the multipath disturbance cannot effectively filter out the white noise of the high frequency part of the signal, and it is prone to edge effect. In this paper, for extracting multipath, a wavelet packet algorithm based on two-dimensional moving weighted average processing (WP-TD) is proposed. This algorithm can not only effectively filter out the white noise of the high frequency part of the signal, but also weaken the influence of the edge effect. Furthermore, considering the repeatability of multipath error in static positioning, we propose a method for determining the level of wavelet packet decomposition layers which make multipath extraction more effectively. The experimental results show that the corrected positioning accuracy is 14.14% higher than that of the traditional wavelet transform when applying the obtained multipath to DD coordinate sequences for position correction.


Introduction
In GPS static high-precision positioning of short baselines, errors such as ionospheric delay, tropospheric delay, satellite orbit error, receiver clock error, and satellite clock error can be eliminated or attenuated by differential techniques while multipath disturbance are not correlated at both ends of the baseline which make it impossible to be eliminated by differential techniques and become the main source of error affecting the positioning accuracy [1,2].
GPS multipath disturbance occurs when GPS signals travel from a satellite to a receiver via several paths due to reflection or diffraction of signals by nearby obstacles [3]. Multipath can distort signal modulation and carrier phase, thereby reducing the accuracy of GPS static high-precision positioning. In addition, this effect may hinder the fixation of ambiguities and even lead to erroneous solutions.
In order to mitigate and eliminate multipath, a variety of different strategies have been proposed, including appropriate site selection, hardware-based methods and software-based methods.
Appropriate site selection is the simplest and the most effective way to solve the multipath effect. The construction of the base station should be chosen in an open area, avoiding various obstacles causing excessive reflection and obtaining as many visible satellite signals as possible.
PozoRuz et al. (1998) proposed a method calculating the positions in the plane using the three highest satellites based on a new satellite selection criterion which can reduce multipath [4].
The hardware-based method refers to improvements of antenna and receiver design.
Scire Scappuzzo et al. (2009) designed a low-multipath wideband GPS antenna with cutoff or non-cutoff corrugated ground plane which can operate uniformly between 1.15 GHz and 1.60 GHz while maintaining the required low-multipath performance in the whole bandwidth [5]. Ray et al. (1999) developed a system for reducing the effect of carrier-phase multipath on static GPS applications using multiple closely spaced antennas [6]. Bryan R. Townsend et al. (1994) introduced a new tracking loop which takes full advantage of the Narrow Correlator spacing design, but, in addition, is much more resistant to multipath effects on the correlation function and thereby reduces the multipath bias on the pseudorange measurements [7]. R.D.J. van Nee et al. (1992) designed a specific receiver structure which simultaneously estimates the parameters of line-of-sight plus multipath signals for reducing GPS code and carrier multipath errors [8].
However, the hardware is inaccessible for most users and it is expensive. Despite this, we cannot deny that hardware methods have better performance in other aspects, such as stability and efficiency. In contrast, software methods are more convenient and cheaper for most people.
The software-based method mainly focuses on the post-processing of signal by using algorithms and filters.
Axelrad et al. (1996) proposed an improved technique which can adaptively estimate the spectral parameters (frequency, amplitude, phase offset) of multipath in the associated signal-to-noise ratio (SNR) [9]. Linlin Ge et al. (2000) proposed an adaptive finite-duration impulse response filter, based on a least-mean-square algorithm, that has been developed to derive a relatively noise-free time series from the CGPS results [10]. There are some other filters, such as Kalman filters (Nce and Sahin 2000) [11], FIR filters (Han and Rizos (1997)) [12], and Vondrak filters (Vondrak (1977) [13], have also been developed to reduce GPS multipath effects.
In addition, wavelet transform has a good localization property within both the time and frequency domains. Particularly, after Stephane Mallat proposed the multi-resolution algorithm, also named Mallat algorithm [14], wavelet transform is widely applied signal processing, including extracting multipath. Huang et al. (2003) applied wavelet transform in dynamic deformation monitoring for high-rise buildings [15]. Zhang and Bartone (2004) applied the wavelet technique to mitigate errors for satellite based navigation systems which can mitigate multipath error in a real-time conditions [16]. E. M. Souza, J. F. G. Monico (2004) apply the wavelet transform to decompose the pseudorange and carrier phase DD signals to separate the high frequencies, which are due to multipath from long delays, and the low frequencies effects, associated with multipath from short delays [17]. P. Zhong et al. (2008) proposed a method based on the technique of cross-validation for automatically identifying wavelet signal layers is developed and used for separating noise from signals in data series, and applied to mitigate GPS multipath effects [3]. M. R. Azarbad and M. R. Mosavi (2014) proposed a new multipath mitigation method based on the wavelet transform (WT). The method uses the stationary wavelet transform (SWT) to decompose the double difference residuals [18]. Lawrence Lau (2017) proposed a generic and robust three-level wavelet packets based denoising method for repeat-time-based carrier phase multipath filtering in relative positioning; the results show that the wavelet packets based method is better than the DWT-based method in the repeat time-based multipath filtering [19]. Souza et al. (2017) proposed a new approach for structure monitoring from GPS multipath effect and wavelet spectrum, and the experience investigated the feasibility of using wavelet spectra analysis of the multipath signal to monitor structure movement [20].
As we know, based on the Mallat algorithm, a signal can be decomposed into scale space and wavelet space. With the increase of the decomposition levels, the scale space would be decomposed continuously, whereas the wavelet space cannot be decomposed, which results in a low resolution in high frequency regions of wavelet transform [21]. This means that the multipath in the high frequency part of the signal cannot be effectively extracted by traditional wavelet transform.
Hence, wavelet packet transform (WPT) was pioneered by Coifman et al. based on wavelet transform [22]. WPT keeps the property of wavelet transform naturally, good localization property within both the time and frequency domains. This property greatly improves the accuracy of WPT in analyzing non-stationary and non-periodic signals. In addition, WPT effectively mitigates the defects in wavelet transform of low resolution in high frequency regions by decomposing both the scale space and wavelet space [21]. This means we can further subdivide the high frequency portion of signal and extract multipath more effectively.
The multipath extraction efficiency partially depends on the chosen mother wavelet for WPT. In this paper, we won't attempt to compare the performance of mother wavelets on retrieving multipath errors from noisy coordinate residual sequences. E. M. Souza et al. (2004) had verified that SYM12, SYM8 and DAUB8 presented the better performance for reconstruction of the GPS DD signal [17]. In this paper, we will choose DAUB8 as the mother wavelet for WPT.
The other challenging issue related to the multipath extraction efficiency is the level of decomposition which depends on the sub-band frequency in which multipath lies [18]. Thence, we also propose a method for determining the number of wavelet packet decomposition levels considering the repeatability of multipath error in static positioning.
In addition, edge effects in signal processing using wavelet packet transform is another un-neglected problem. P. Zhong et al. (2008) choose about 70% of the data in the middle of the observational series for cross-validation to prevent edge effects due to poorer filtering results at the ends of a data series [3]. This method solves the edge effect problem in a certain condition, but, when the coordinate sequence residual is short, it will become invalid. Therefore, referring to the idea of moving weighted average method and the principle of bilateral filters in image processing [23,24], we propose a new algorithm named the two-dimensional moving weighted average algorithm. The double-difference coordinate residual sequences are smoothed by the proposed algorithm to achieve the purpose of edge-preserving and pre-denoising.
In general, we propose a new method named wavelet packet algorithm based on two-dimensional moving weighted average processing (TDMWA), compared with the traditional wavelet algorithm, which can not only more effectively mitigate the multipath of the DD coordinate residual sequences, but also effectively weaken the influence of the edge effect.
This study is organized as follows. The principle of position correction based on multipath periodicity is introduced in the Section 2. The theory of the two-dimensional moving weighted average algorithm and the preprocessing of smoothing the double-difference coordinate sequence residual obtained by static observation for three consecutive days using TDWMA algorithm are introduced in Section 3. The theory of wavelet packet transform, the theory of correlation analysis of smoothed double-difference coordinate residual sequences, the specific steps for determining the number of wavelet packet decomposition layers and the principle of wavelet packet threshold denoising will be introduced in Section 4. The overall program of mitigating the multipath in GPS static high-precision positioning will be introduced in Section 5. The experimental verification will be introduced in Sections 6 and 7. The data set used in Section 6 is the simulation data for verifying the feasibility of the proposed method. Section 7 uses the GPS data acquisition system built to collect the GPS data to further verify the feasibility of the actual environment and the denoising performance of the Vondrak algorithm, traditional WT algorithm and the WP-TD algorithm is compared from the correlation and root mean square error. Conclusions are given in Section 8.

The Principle of Position Correction Based on Multipath Periodicity
In this paper, we take the double difference carrier phase observation residuals obtained from the reference day (Day 1) as the correction to the observations on the adjacent days (Day 2 and Day 3) considering that multipath is repetitive in static observations of the fixed point. Methods for position correction based on multipath repeatability can be found in Khelifa et al. (2011) [25], Ye et al. (2013) [26], and Lawrence Lau et al. (2017) [19]. The double difference carrier phase observations [27] indicated as Equation (1): where Φ nm rb denotes the double difference carrier phase observation between satellites n and m, and stations r and b. ρ nm rb denotes the geometric distance from the satellite center to the antenna phase center. N nm rb denotes the double difference integer ambiguity. ε nm φ,rb denotes the double difference carrier phase measurement noise. M nm φ,rb denotes the double difference carrier phase multipath. λ denotes the carrier wavelength.
In GPS high-precision relative positioning of short baseline (less than 3 km). Since the DD carrier phase multipath error is always less than a quarter of the carrier wavelength, the observation residuals obtained on the reference day do not take the double difference integer ambiguity into consideration [21].
Based on the double difference carrier phase observations, the positioning solution of each epoch on the reference day indicated as Equation (2): where .
X denotes the best estimated positioning solution of the fixed point in static observation. X denotes the known positioning solution of the fixed point in static observation. M denotes the multipath, and ε denotes the observation noise (white noise). Subscript Re indicates the reference day.
Moving the known item of Equation (2) to the left side to get Equation (3): The right side of the Equation (3) is the positioning residuals that consist of multipath and observation noise (white noise), which will be regarded as the DD coordinate residual sequences in the next section.
Because the repeat-time-based multipath increases the observation noise level in the positioning solutions of the fixed point on the adjacent days, it is necessary to eliminate the noise of the DD coordinate residual sequences and obtain the noiseless multipath as the correction of the positioning solutions of the adjacent days.
The corrected position of each epoch on the adjacent day indicated as Equation (4): whereX denotes the multipath-corrected positioning solutions. Subscript Ad indicates the adjacent days.

The Basic Principle of the TDMWA Algorithm
In GPS static high-precision positioning, the multipath tends to be stable between the adjacent epochs, but the DD coordinate residual sequences fluctuation intensified between adjacent epochs due to the existence of white noise. Even worse, there exist poorer filtering results at the end of the DD coordinate residual sequences because of the edge effects.
In order to weaken the influence of edge effects and weaken the white noise, we propose a new method named two-dimensional moving weighting average algorithm (TDMWA) for smoothing the DD coordinate residual sequences. Considering that the white noise is independent and obeys the normal distribution, the algorithm performs moving weighted average processing on each of the  α j is the weighting coefficient of the j-th delay epoch in the moving average period T, and β k is the weighting coefficient of the k-th decentralized epoch in the moving average period N r , then where x i represents the all observation epoch coordinates, and i ∈ I. Then, calculating the standard deviation of the observation epoch coordinate, and using the standard deviation as the threshold θ I for selecting decentralized epoch coordinate constituting Q: Then, the value of the current observation epoch coordinate x r becomes the new mean of all the observation epochs, Next, calculating the standard deviation σ i of all the observation epochs coordinate x i according to the new mean obtained in Equation (9), The observation epoch coordinate whose standard deviation σ i is not greater than the threshold θ I is used as the decentralized epoch coordinate x k , which will be selected to determine the smoothed value of the r-th observation epoch in the value domain, and the number of the decentralized epoch coordinate x k is the moving average period N r , and N r ≤ I.
In the value domain, the difference between each decentralized epoch coordinate x k and the current observation epoch coordinate x r is calculated, and the weighting coefficient are assigned to the respective decentralized epoch coordinate x k according to the magnitude of the difference, and, the smaller the difference, the greater the influence of similarity effect, the greater the weight of the corresponding decentralized epoch coordinate x k . Let µ be the average of the magnitude of the difference of the decentralized epoch coordinate x k , ∆µ k is the difference of the each decentralized epoch coordinate x k with respect to the current observation epoch coordinate x r , and ϑ k is the variation coefficient of the each decentralized epoch coordinate x k . Then, the weighting coefficient of the decentralized epoch coordinate x k in the current epoch r can be calculated by the following formula: Then, we substitute T, N r , x j , x k , α j , β k into Equation (5) where Υ = X T r +X Nr Finally, we substitute the obtained γ T and γ N r into Equation (13) to obtain the new valueX r of the current observation epoch r, which is the smoothed observation epoch coordinate.
Following the steps as above, this method is expected to solve the edge effect and achieve the purpose of edge preservation and denoising, but using this method will retain too much high frequency information, so that the white noise in high frequency will not be completely eliminated, so this paper will add wavelet packet transform based on this method to further eliminate the white noise in high frequency.

The Theory of WPT
From the introduction, we have known that wavelet packet transform provides a more sophisticated analysis of the signal. Wavelet packet transform divides the time-frequency space into more detail, and it has higher resolution for the high-frequency part of the signal than the binary wavelet transform. Moreover, based on the theory of wavelet analysis, it introduces the concept of optimal basis selection. That is, after the frequency band is divided into multiple levels, according to the characteristics of the analyzed signal, the optimal basis function is adaptively selected to match the signal to improve the signal analysis capability. In order to further eliminate the white noise of the high frequency part of the smoothed DD coordinate residual sequences, which wavelet packet transform is applied to decompose, the wavelet packet decomposition algorithm [21] can be indicated as Equation (15): where d 2ω−1 M denotes the 2ω − 1 wavelet packet coefficient of the M-th layer, d ω M−1 represents the ω-th wavelet packet coefficient on the M-1 layer, and p 2ω M represents the 2ω-th wavelet packet coefficients on the M-th layer, G is the decomposition filter relation to scale function, H is the decomposition filter of wavelet function, and M ≥ 1, The wavelet packet reconstruction algorithm [21] can be indicated as Equation (16): where d 2ω−1 M+1 represents the 2ω − 1 wavelet packet coefficient on the M + 1 layer, d 2ω M+1 indicates the 2ω wavelet packet coefficients on the J + 1 layer, g is the reconstruction filter related to scaling function, and h is the reconstruction filter related to wavelet function.
From the perspective of time-frequency conversion, the DD coordinate residual sequences after the first WPT decomposition could be converted into two equal frequency bands by Equation (15) continuously. An example of three level decomposition of WPT is shown as Figure 1. can be divided continuously.
An example of three level decomposition of WPT is shown as Figure 1. In Figure 1, A represents a low frequency, D represents a high frequency, and the serial number at the end represents the number of layers of wavelet decomposition.
The decomposition relationship can be expressed as Equation (17): After the DD coordinate residual sequences are decomposed according to the structure shown in Figure 1, using the thresholding denoising mentioned in the Section 4.4, the wavelet packet coefficients that satisfy the condition will be used to reconstruct the compressed DD coordinate residual sequences while others will be set to zero.

The Theory of Correlation Analysis
In view of the strong correlation of multipath in the same time between the adjacent days in static positioning, though the smoothed DD coordinate residual sequences still contain residual white noise, but the multipath error is dominant, which contributes to strong correlation between the smoothed DD coordinate residual sequences.
We have the decomposed DD coordinate residual sequences in Section 3.1, which can be applied to correlation analysis for determining the best decomposition level of the WPT.
In the three consecutive days, regard one day as the reference day and the rest as the adjacent days.
A is the wavelet packet coefficient of the low-frequency part of the decomposed DD coordinate In Figure 1, A represents a low frequency, D represents a high frequency, and the serial number at the end represents the number of layers of wavelet decomposition.
The decomposition relationship can be expressed as Equation (17): After the DD coordinate residual sequences are decomposed according to the structure shown in Figure 1, using the thresholding denoising mentioned in the Section 4.4, the wavelet packet coefficients that satisfy the condition will be used to reconstruct the compressed DD coordinate residual sequences while others will be set to zero.

The Theory of Correlation Analysis
In view of the strong correlation of multipath in the same time between the adjacent days in static positioning, though the smoothed DD coordinate residual sequences still contain residual white noise, but the multipath error is dominant, which contributes to strong correlation between the smoothed DD coordinate residual sequences.
We have the decomposed DD coordinate residual sequences in Section 3.1, which can be applied to correlation analysis for determining the best decomposition level of the WPT.
In the three consecutive days, regard one day as the reference day and the rest as the adjacent days.
A is the wavelet packet coefficient of the low-frequency part of the decomposed DD coordinate sequences of the reference day, and B is the wavelet packet coefficient of the corresponding low-frequency part of the decomposed DD coordinate sequence of the adjacent day, and ρ AB is the cross-correlation coefficient.
The formula for calculating the correlation coefficient can be expressed as Equation (18): where The magnitude of ρ AB reflects the correlation degree between the analyzed sequences. The closer the correlation coefficient ρ AB is to 1, the closer A and B are to linear correlation, and the greater the correlation of A and B.
However, due to the influence of residual white noise, satellite visible conditions and the distance between the antenna of the monitoring station and the reflecting surface, the correlation coefficient of the smoothed DD coordinate residual sequences between adjacent days becomes smaller.

The Method for Determining the Best Decomposition Level of the WPT
Determining the best decomposition level of the WPT that is related to the validity of multipath extraction, if there are too many decomposition layers, the useful information in the smoothed DD coordinate residual sequences will be lost. If the number of decomposition layers is too small, the white noise in the smoothed DD coordinate residual sequences cannot be eliminated absolutely.
The method based on the correlation analysis for determining the best decomposition level of the WPT can be described as follows: Step1: Calculating the cross-correlation coefficient of the smoothed DD coordinate residual sequences between the reference day and other adjacent days, for any reference day, and selecting the smoothed DD coordinate residual sequences that have the largest cross correlation coefficients.
Step2: Performing M-layers decomposition on the set of the smoothed DD coordinate residual sequences obtained in step 1 by using the WPT algorithm, and calculating the cross-correlation coefficient of the wavelet packet coefficients of the approximate part (low-frequency part) of the each decomposition layer between the smoothed DD coordinate residual sequences obtained in step 1.
Step3: Determining whether the target layer J is smaller than M, where M is the preset number of decomposition layers, and the target layer J is the number of layers in which the maximum value of the cross-correlation coefficient is located.
If yes, determining the target layer J be the best decomposition level of the WPT; If not, the value of M is incremented by one, and the process returns to step 3 until the target layer J is smaller than M. By this way, we can figure out the best decomposition level of the WPT that can be applied to analyze the smoothed DD coordinate residual sequences in the optimal state.

The Principle of Wavelet Packet Threshold Denoising
We have the decomposed DD coordinate residual sequences with the best decomposition level of the WPT and the wavelet packet coefficients obtained should be denoised by thresholding for reconstructing the DD coordinate residual sequences.
The process of denoising by thresholding is carried out by comparing the magnitude of the wavelet packet coefficients d i with a threshold λ. The process can be described in two aspect: one is the choice of the threshold function F and the other is the choice of the threshold parameter λ.
There are three choices of threshold function F including hard thresholding F H λ (d i ), soft thresholding F S λ (d i ) and quantitative thresholding F Q λ (d i ) [28]: where d i denotes the wavelet coefficients, and P is the value for which a certain percentage of coefficients d i is eliminated [17]. In this paper, the hard thresholding function is chosen for denoising the wavelet packet coefficients, which had been proved be the most suitable function for GPS applications [17].
As for the choice of the threshold parameter λ, Donoho (1995) had proposed the universal threshold [29]: where n denotes the number of samples in white noise vector N and σ is the white noise level and needs to be estimated in each DD coordinate residual sequence. For determining the value of σ, Donoho and Johnstone (1994) proposed the following estimator [30]: The noise vector N in the first level of decomposition layers is used for determining the value of σ.
Applying the thresholding method mentioned as above, all of the wavelet packet coefficients are denoised, which contribute to the reconstruction of the DD coordinate residual sequences and obtain the noiseless multipath.

The Overall Program
The overall program of mitigating the multipath in GPS static high-precision positioning is shown in Figure 2. denoised, which contribute to the reconstruction of the DD coordinate residual sequences and obtain the noiseless multipath.

The Overall Program
The overall program of mitigating the multipath in GPS static high-precision positioning is shown in Figure 2.

Verify the Feasibility of the Algorithm under Simulated Conditions
The DD coordinate residual sequences is simulated with the following model: where ( ) consists of a cosine signal with a period of 1200 s and two sinusoidal signals with periods of 900 s and 300 s which represents typical GPS multipath wavelengths, and ( ) denotes a Gaussian white noise sequence. Take the noise level of ( ) as N (0, 1.5 2 ). The data sampling rate is 1 s, and the sample size is 6000. The simulated signal is denoised by WT transform, the TDMWA algorithm and WP-TD algorithm, respectively. DAUB8 is taken as the wavelet function, and we assume that the best decomposition level is 5.
The denoising result is shown in Figure 3, and the relative distribution of the signal denoised by the three different algorithms with the original signal is shown in Figure 4.

Verify the Feasibility of the Algorithm under Simulated Conditions
The DD coordinate residual sequences is simulated with the following model: where S(t) consists of a cosine signal with a period of 1200 s and two sinusoidal signals with periods of 900 s and 300 s which represents typical GPS multipath wavelengths, and e(t) denotes a Gaussian white noise sequence. Take the noise level of e(t) as N (0, 1.5 2 ). The data sampling rate is 1 s, and the sample size is 6000. The simulated signal is denoised by WT transform, the TDMWA algorithm and WP-TD algorithm, respectively. DAUB8 is taken as the wavelet function, and we assume that the best decomposition level is 5.
The denoising result is shown in Figure 3, and the relative distribution of the signal denoised by the three different algorithms with the original signal is shown in Figure 4.
The residual signal generated by the three different algorithms is shown in Figure 5. From Figures 4 and 5, we can visually see that the signal denoised by the WP-TD algorithm better restores the original signal. The denoised result of the signal processed using the TDMWA algorithm is also slightly better than the traditional WT algorithm. In addition, the signal denoised using the WT algorithm has an edge effect, but the signal denoised using the TDMWA algorithm does not appear.
The signals denoised by the above three algorithms and the residual signals generated are analyzed quantitatively. The analysis indicators include the root mean square (RMS) value of the noise part of denoised signal N RMS , the RMS value of the signal part of the denoised signal S RMS and the correlation coefficient ρ between the denoised signal and the original signal. In order to scientifically evaluate the performance of the three algorithms, simulation analysis was carried out in four different Gaussian white noise simulation environments. The results are shown in Tables 1-3.  The residual signal generated by the three different algorithms is shown in Figure 5.  The residual signal generated by the three different algorithms is shown in Figure 5. From Figures 4 and 5, we can visually see that the signal denoised by the WP-TD algorithm better restores the original signal. The denoised result of the signal processed using the TDMWA algorithm is also slightly better than the traditional WT algorithm. In addition, the signal denoised using the WT algorithm has an edge effect, but the signal denoised using the TDMWA algorithm does not appear.
The signals denoised by the above three algorithms and the residual signals generated are analyzed quantitatively. The analysis indicators include the root mean square (RMS) value of the noise part of denoised signal , the RMS value of the signal part of the denoised signal and the correlation coefficient between the denoised signal and the original signal. In order to scientifically evaluate the performance of the three algorithms, simulation analysis was carried out in four different Gaussian white noise simulation environments. The results are shown in Tables 1-3.   On one hand, it can be seen from Tables 1-3 that the correlation coefficients of the three algorithms mostly exceed 0.95 at different noise levels, indicating that all of them can restore the information of the original very well. However, the value of S RMS of the WP-TD algorithm is relatively optimal at different noise levels, and its correlation coefficient is greater than 0.98.
On the other hand, the value of N RMS of the three algorithms are basically equivalent while there is some difference in the value of S RMS . According to the tables, as the noise level increases, the value of S RMS of the WT algorithm and the TDMWA algorithm increase significantly compared with WP-TD algorithm. However, at high noise levels, the TDMWA algorithm increases slowly while the WT algorithm is relatively fast, and its S RMS is still superior to the WT algorithm.
In summary, we can draw such a conclusion that the WP-TD algorithm can denoise the noise signal more effectively, and can effectively weaken the influence of the edge effect in signal filtering.

Analyze the Performance of the Algorithm in a Measured Environment
In order to further verify the ability of the WP-TD algorithm to eliminate noise in the actual environment, we established a GPS data acquisition system in the laboratory and conducted experiments in a multipath environment.

Set up GPS Data Acquisition System
The block diagram of the overall scheme and key modules of the positioning data acquisition system are shown in Figures 6 and 7, which shows the scene map and physical device diagram.
In this system, the base station and mobile station are all equipped with a DC-RTK-00 Beidou/GPS/Galileo three-mode single-frequency RTK module, a power model (5 V), a microcontroller (Raspberry3 b+) and a data transmission module (Ethernet card/4G card). The base station sends the differential data to the server over the wired network and is cached by the lab server. The mobile station accesses the server through the 4G network to obtain differential data, and, at the same time, the mobile station equipped with the 4G network card transmits the positioning result back to the lab server.  On one hand, it can be seen from Tables 1-3 that the correlation coefficients of the three algorithms mostly exceed 0.95 at different noise levels, indicating that all of them can restore the information of the original very well. However, the value of of the WP-TD algorithm is relatively optimal at different noise levels, and its correlation coefficient is greater than 0.98.
On the other hand, the value of of the three algorithms are basically equivalent while there is some difference in the value of . According to the tables, as the noise level increases, the value of of the WT algorithm and the TDMWA algorithm increase significantly compared with WP-TD algorithm. However, at high noise levels, the TDMWA algorithm increases slowly while the WT algorithm is relatively fast, and its is still superior to the WT algorithm. In summary, we can draw such a conclusion that the WP-TD algorithm can denoise the noise signal more effectively, and can effectively weaken the influence of the edge effect in signal filtering.

Analyze the Performance of the Algorithm in a Measured Environment
In order to further verify the ability of the WP-TD algorithm to eliminate noise in the actual environment, we established a GPS data acquisition system in the laboratory and conducted experiments in a multipath environment.

Set up GPS Data Acquisition System
The block diagram of the overall scheme and key modules of the positioning data acquisition system are shown in Figures 6 and 7, which shows the scene map and physical device diagram.
(a) Overall plan

Data Collection and Processing
The fixed point is located at the top of the research building of Beijing University of Posts and Telecommunications. The settings of related parameters of data acquisition are shown in Table 4. The DD coordinate residual sequences to be processed are part of the elevation data sequence measured of the fixed point for three consecutive days in a multipath environment. Figure 8 shows the DD coordinate residual sequences of the elevation direction for three consecutive days. As can be seen from Figure 8, under the influence of the observation noise, the DD coordinate residual sequences for three consecutive days still have significant repeatability.
Beidou/GPS/Galileo three-mode single-frequency RTK module, a power model (5 V), a microcontroller (Raspberry3 b+) and a data transmission module (Ethernet card/4G card). The base station sends the differential data to the server over the wired network and is cached by the lab server. The mobile station accesses the server through the 4G network to obtain differential data, and, at the same time, the mobile station equipped with the 4G network card transmits the positioning result back to the lab server.

Data Collection and Processing
The fixed point is located at the top of the research building of Beijing University of Posts and Telecommunications. The settings of related parameters of data acquisition are shown in Table 4. The DD coordinate residual sequences to be processed are part of the elevation data sequence measured of the fixed point for three consecutive days in a multipath environment. Figure 8 shows the DD coordinate residual sequences of the elevation direction for three consecutive days. As can be seen from Figure 8, under the influence of the observation noise, the DD coordinate residual sequences for three consecutive days still have significant repeatability.  According to the overall program of mitigating the multipath in GPS static high-precision positioning established in this paper, we need to smooth the DD coordinate residual sequences first using the TDMWA algorithm.
Secondly, we need to analyze the correlation of the DD coordinate sequences residuals for three consecutive days and determine the best decomposition level of WPT based on the method proposed in this paper.
The process of determining the best decomposition level of WPT is shown in Figure 9.
According to the overall program of mitigating the multipath in GPS static high-precision positioning established in this paper, we need to smooth the DD coordinate residual sequences first using the TDMWA algorithm.
Secondly, we need to analyze the correlation of the DD coordinate sequences residuals for three consecutive days and determine the best decomposition level of WPT based on the method proposed in this paper.
The process of determining the best decomposition level of WPT is shown in Figure 9. From right to left, and from top to bottom are the DD coordinate residual sequence of the reference day preprocessed by the TDMWA algorithm, the seven-layer decomposition tree structure of the DD coordinate residual sequence of the reference day, the cross-correlation colored coefficients of the corresponding nodes and the value of the first wavelet packet coefficient of the decomposed DD coordinate residual sequence of the reference day located on the sixth layer. It can be seen from the distribution of the colored coefficients in the figure that the larger values of the correlation coefficient are mostly distributed in the sixth layer and the fifth layer, and the maximum value appears in the sixth layer.
According to the judgment conditions in the method proposed in this paper, we can determine that the best denoise result can be obtained when the decomposition level of the wavelet packet is six layers.
Next, we perform a six-layer wavelet packet decomposition on the smoothed DD coordinate residual sequences and reconstruct the decomposed DD coordinate residual sequences based on the wavelet packet coefficients obtained by threshold denoising.
Finally, we can obtain almost noiseless multipath, which can be applied to correct the position sequences of the adjacent days.
The DD coordinate residual sequences denoised by the WP-TD algorithm, Vondrak algorithm and WT algorithm for three consecutive days are shown in Figure 10, and Figure 11 shows the difference of denoised DD coordinate residual sequences for the second day and the third day compared to the first day, respectively, and Figure 12 shows the noise residuals denoised by the three algorithms on the first day. From right to left, and from top to bottom are the DD coordinate residual sequence of the reference day preprocessed by the TDMWA algorithm, the seven-layer decomposition tree structure of the DD coordinate residual sequence of the reference day, the cross-correlation colored coefficients of the corresponding nodes and the value of the first wavelet packet coefficient of the decomposed DD coordinate residual sequence of the reference day located on the sixth layer. It can be seen from the distribution of the colored coefficients in the figure that the larger values of the correlation coefficient are mostly distributed in the sixth layer and the fifth layer, and the maximum value appears in the sixth layer.
According to the judgment conditions in the method proposed in this paper, we can determine that the best denoise result can be obtained when the decomposition level of the wavelet packet is six layers.
Next, we perform a six-layer wavelet packet decomposition on the smoothed DD coordinate residual sequences and reconstruct the decomposed DD coordinate residual sequences based on the wavelet packet coefficients obtained by threshold denoising.
Finally, we can obtain almost noiseless multipath, which can be applied to correct the position sequences of the adjacent days.
The DD coordinate residual sequences denoised by the WP-TD algorithm, Vondrak algorithm and WT algorithm for three consecutive days are shown in Figure 10, and Figure 11 shows the difference of denoised DD coordinate residual sequences for the second day and the third day compared to the first day, respectively, and Figure 12 shows the noise residuals denoised by the three algorithms on the first day. Sensors 2019, 19, x FOR PEER REVIEW 17 of 21 Figure 10. The DD coordinate residual sequences denoised by the three algorithms for three consecutive days. Figure 11. The difference of DD coordinate residual sequences denoised by the three algorithms.   It can be seen from the Figures 10 and 11 that the difference of DD coordinate residual sequence processed by the WP-TD algorithm is much smaller than the WT algorithm and the Vondrak algorithm, which shows that applying the DD coordinate residual sequence denoised by the WP-TD algorithm to correct the position sequence of the adjacent day will be more precise.
Moreover, the DD coordinate residual sequence has a larger extremum at the end of sequence, which was denoised by the WT algorithm without the preprocessing of TDMWA algorithm, while the DD coordinate residual sequence presents a smooth curve at the end of the sequence which is preprocessed by the TDMWA algorithm before being denoised by the WPT algorithm.

Data Analysis
In order to further compare the performance for extracting multipath of the three algorithms in the actual environment, the RMS value of the DD coordinate residual sequences before and after filtering, the RMS value of the DD coordinate residual sequences of the adjacent days corrected by the almost noiseless multipath obtained in the reference day and the correlation coefficient of the DD coordinate residual sequence of the three days are used as the evaluation indicator.
The comparison results of the three algorithms are shown in Tables 5-7, respectively. First of all, it can be seen from the results in Table 5 that the RMS value of the DD coordinate residual sequence after filtering by the WP-TD algorithm is basically equivalent to the RMS value of the DD coordinate residual sequence after filtering by the Vondrak algorithm, but both are slightly better than the RMS value of the DD coordinate residual sequence after filtering by the WT algorithm. In addition, the noise residuals denoised by the three algorithms on the first day are not much different as shown in Figure 12, indicating that the observation noise (white noise) is only a small part of the DD coordinate residuals sequence, and multipath dominates. It can be seen from the Figures 10 and 11 that the difference of DD coordinate residual sequence processed by the WP-TD algorithm is much smaller than the WT algorithm and the Vondrak algorithm, which shows that applying the DD coordinate residual sequence denoised by the WP-TD algorithm to correct the position sequence of the adjacent day will be more precise.
Moreover, the DD coordinate residual sequence has a larger extremum at the end of sequence, which was denoised by the WT algorithm without the preprocessing of TDMWA algorithm, while the DD coordinate residual sequence presents a smooth curve at the end of the sequence which is preprocessed by the TDMWA algorithm before being denoised by the WPT algorithm.

Data Analysis
In order to further compare the performance for extracting multipath of the three algorithms in the actual environment, the RMS value of the DD coordinate residual sequences before and after filtering, the RMS value of the DD coordinate residual sequences of the adjacent days corrected by the almost noiseless multipath obtained in the reference day and the correlation coefficient of the DD coordinate residual sequence of the three days are used as the evaluation indicator.
The comparison results of the three algorithms are shown in Tables 5-7, respectively.  First of all, it can be seen from the results in Table 5 that the RMS value of the DD coordinate residual sequence after filtering by the WP-TD algorithm is basically equivalent to the RMS value of the DD coordinate residual sequence after filtering by the Vondrak algorithm, but both are slightly better than the RMS value of the DD coordinate residual sequence after filtering by the WT algorithm. In addition, the noise residuals denoised by the three algorithms on the first day are not much different as shown in Figure 12, indicating that the observation noise (white noise) is only a small part of the DD coordinate residuals sequence, and multipath dominates.
Secondly, the results in Table 6 show that the GPS multipath of the fixed point is highly repetitive, and the correlation coefficient increases after filtering observation noise. Moreover, the performance of filtering of the WP-TD algorithm is better than WT algorithms and Vondrak algorithms more or less. Through calculation, we can know that the correlation of the DD coordinate residual sequence of the three days filtered by WP-TD is 3.02% and 1.78% higher than that of WT algorithms and Vondrak algorithms, respectively.
Last but not least, according to the data in Table 7, after correcting by the denoised multipath of the reference day (Day 1) which was extracted by the WT algorithms, the Vondrak algorithms, and the WP-TD algorithm, we can figure out that the accuracy of the position sequence of the adjacent days (Day 2 and Day 3) is increased by 57.47%, 65.98%, and 71.61%, respectively, compared with before correction, which further verified that the WP-TD algorithm can obtain an more accurate multipath correction model.

Conclusions
The WP-TD algorithm is proposed, in which the method of determining the optimal decomposition level of WPT is given, which can avoid the redundancy and loss of information after signal reconstruction to the greatest extent.
The WP-TD algorithm is proposed, which makes use of the TDMWA algorithm, and is proven to have good performance that weakens the influence of the edge effect through simulation analysis and measured data analysis.
The WP-TD algorithm is proposed, which uses the WPT to divide the signal frequency band into multiple layers, further decomposes the high frequency part without subdivision in the traditional WT algorithm, and eliminates the observation noise in the high frequency part.
The WP-TD algorithm is proposed, which applies the denoised multipath of the reference day to correct the position sequence of the adjacent days, improving positioning accuracy by 14.4%, compared to the traditional WT algorithm.
In summary, the WP-TD algorithm is proposed, which can perform well for eliminating the multipath in GPS static high-precision positioning.