1. 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
Section 6 and
Section 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.
2. 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
denotes the double difference carrier phase observation between satellites
n and
m, and stations
r and
b.
denotes the geometric distance from the satellite center to the antenna phase center.
denotes the double difference integer ambiguity.
denotes the double difference carrier phase measurement noise.
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
denotes the best estimated positioning solution of the fixed point in static observation.
denotes the known positioning solution of the fixed point in static observation.
denotes the multipath, and
denotes the observation noise (white noise). Subscript
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):
where
denotes the multipath-corrected positioning solutions. Subscript
indicates the adjacent days.
4. Time-Frequency Conversion by WPT
4.1. 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
denotes the 2ω − 1 wavelet packet coefficient of the
M-th layer,
represents the
ω-th wavelet packet coefficient on the M-1 layer, and
represents the 2
ω-th wavelet packet coefficients on the
M-th layer,
is the decomposition filter relation to scale function,
is the decomposition filter of wavelet function, and
.
The wavelet packet reconstruction algorithm [
21] can be indicated as Equation (16):
where
represents the
wavelet packet coefficient on the
layer,
indicates the
wavelet packet coefficients on the
layer,
is the reconstruction filter related to scaling function, and
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), i.e., and , where is Nyquist frequency. The WPT coefficients are divided into and as well. Moreover, both and 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.
4.2. 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.
is the wavelet packet coefficient of the low-frequency part of the decomposed DD coordinate sequences of the reference day, and is the wavelet packet coefficient of the corresponding low-frequency part of the decomposed DD coordinate sequence of the adjacent day, and is the cross-correlation coefficient.
The formula for calculating the correlation coefficient can be expressed as Equation (18):
where
is the covariance of variables
and
,
is the variance of the variable
,
is the variance of the variable
,
is the average value of the variable
,
is the average value of the variable
.
The magnitude of reflects the correlation degree between the analyzed sequences. The closer the correlation coefficient is to 1, the closer A and B are to linear correlation, and the greater the correlation of and .
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.
4.3. 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.
4.4. 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 with a threshold . The process can be described in two aspect: one is the choice of the threshold function and the other is the choice of the threshold parameter .
There are three choices of threshold function
including hard thresholding
, soft thresholding
and quantitative thresholding
[
28]:
where
denotes the wavelet coefficients, and
is the value for which a certain percentage of coefficients
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
denotes the number of samples in white noise vector
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 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.
6. 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.52). 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
Figure 4 and
Figure 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
Table 1,
Table 2 and
Table 3.
On one hand, it can be seen from
Table 1,
Table 2 and
Table 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.
7. 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.
7.1. 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
Figure 6 and
Figure 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.
7.2. 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.
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.
It can be seen from the
Figure 10 and
Figure 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.
7.3. 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.
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.
8. 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.