Abstract
Although GPS kinematic relative positioning can provide high accuracy, GPS observables, like any other kind of measurement, are not free of errors. Indeed, they have several kinds of errors. In this paper, we show how to construct a functional mathematical model within the context of a Kalman Filter in order to eliminate most of these errors. Furthermore, we discuss how the multipath effect, a kind of error not modeled in the functional model, can be corrected using the proposed wavelet method. The behavior of the double difference functional model in the kinematic mode is also demonstrated and analyzed aiming to provide better insight into the problem. The results obtained from the multipath experiments were very promising and are presented here.
1. Introduction
Global Navigation Satellite Systems (GNSSs), in particular Global Positioning Systems (GPSs), are revolutionizing navigation as a basic technique for positioning. In GPS relative positioning, the objective is to determine the coordinates of an unknown point with respect to other known point(s). This means that a vector between the two points, which is called the baseline vector or simply the baseline, is determined. To perform this type of positioning, two kinds of GPS observables can be used: pseudoranges (PR) and carrier phases (CP). Like other survey measurements, GPS observables are not free of errors. The increasing demand for high accuracy has required an understanding of the error sources involved as well as of the methods that may reduce or eliminate these errors. Among the error sources, the more significant ones are those from receiver clocks, satellite clocks, the ionosphere, the troposphere, the orbit, and the multipath [1–3].
In relative positioning, single and double differences (DDs) of GPS observables are commonly used to construct the functional model, as they can eliminate or reduce most GPS errors. This functional model reduces most errors if the baseline is short, that is, no more than 20 km. However, the multipath error is the only one that is not eliminated even for short baselines. This is because the multipath depends on the geometry and environment of each point where the GPS antenna is collecting the signals (observables). Therefore, multipath is a major residual error source in double-differenced GPS observables. In kinematic positioning, the nonstationary behavior of the multipath effect is worse, and it is very difficult to remove it from the data.
Although many studies have attempted to mitigate multipath for kinematic positioning, this effect remains a challenge for research.
Multipath effects and signal blockage in GPS navigation in the vicinity of the International Space Station (ISS) were analyzed by [4]. The results showed that the combined effects of GPS blockage and the multipath make GPS-only navigation at the ISS very difficult, if not impossible.
Receiver [5, 6] and hardware technologies [7] have been developed to reduce the multipath, but this effect is only partially mitigated.
Multiple proposals have been put forward with regard to processing techniques, such as those using code-minus-phase measurements that accurately separate or eliminate the multipath signals [8]. However, methods to reduce multipath signals using CP measurements are rare. Reference [9] developed a Receiver Autonomous Integrity Monitoring (RAIM) algorithm by applying dual L1/L2 frequencies to exclude multipath satellites. The algorithm detects multipath satellites by comparing CP with code-minus-phase measurements. However, in this case it is necessary to have a receiver that collects the GPS observables in two bands of frequencies, L1 and L2. A technique that can be applied to either L1 or L1/L2 observables would be very useful, mainly because L1 receivers are the most common due to their lower price.
Spectral analysis has a powerful technique to analyze this kind of nonstationary signal: the wavelet transform [10]. Even so, some processes and specific processing methods must be used together in order to detect and efficiently mitigate the multipath. The signal is analyzed in small windows of data during kinematic processing. To process data for each instant, data from some previous instants are used in multipath detection, but only the observations of the current instant are corrected. In this paper, experiments were carried out in a kinematic mode with a controlled and known vehicle movement. The data were collected in the presence of a reflector surface placed close to the vehicle to cause a multipath. The DD residuals in each instant were analyzed and compared. The obtained results show that the proposed method is promising and efficient in detecting and correcting the multipath effect from both PR and CP DD observables, and consequently, in improving the functional mathematical model.
In summary, we show how errors are accounted for in the functional model construction in relative positioning as well as how the remaining multipath error is detected and corrected by the proposed method using wavelet techniques. Furthermore, the functional models are described, and their behaviors in the kinematic mode are demonstrated and analyzed.
2. Relative Positioning Models
To perform GPS relative positioning, the Kalman Filter (KF) is usually applied. KF deals with two important components: a mathematical functional model and a stochastic model. The functional model describes the mathematical relationship between the GPS observables and the unknown parameters, while the stochastic model describes the statistical characteristics of the GPS observables. The latter therefore depends on the choice of the functional model. The DD technique is commonly used to construct the functional model. In stochastic models, it is usually assumed that all measurements of each observable (PR or CP) have equal variance, and that they are statistically independent [11, 12]. However, stochastic models weighted by satellite elevation angles are also used.
In the next subsections, models of PR and CP observables are described, together with how the errors involved in these observables are eliminated in the functional model.
2.1. GPS Observables
There are two important GPS observables: PR and CP. Measuring PRs and CPs involves advanced techniques in electronics and signal processing [2].
PR is related to the distance between the satellite and the receiver’s antenna, implied by the instants of emission and reception of the PRN codes. PR measurements are obtained from correlation of the code generated by the satellite in the transmission instant (tt) with those generated internally by the receiver at the reception instant [1, 13]. Furthermore, the mathematical expression for the PR between a satellite s and a receiver r must take into account the errors involved:
where
(i) is the Euclidian distance, also called geometric distance, in meters, between the satellite s in tt and the receiver r in tr: with the vectors and representing the Cartesian coordinates of the satellite and the receiver, respectively;(ii) is the error in , in meters, due to errors in the satellite orbit;(iii) and dtr are the satellite and receiver clock errors, in seconds, in relation to the GPS time in the transmission and the reception time, respectively;(iv) is the velocity of light in a vacuum, in m/s;(v) and are the effects, in meters, due to ionospheric and tropospheric influences on the GPS signal;(vi) is the multipath error for PR, in meters; (vii) is the PR measurement noise and unmodelled errors, between 0.5 and 1 meters, depending on the technology used.The CP measurement at a nominal frequency f, in cycles, is determined from the difference between the phase generated by the satellite (s) in tt and its replica generated by the receiver (r) in tr [12, 13]. The CP can be scaled to a unit of length in meters by multiplying it by :
where Nrs is the integer ambiguity. This integer refers to the first measurements and remains constant during the whole period of observation as long as there is no signal loss. During this period, the receiver accumulates the phase differences between arriving phases and internally generated ones. Therefore, the receiver keeps an accumulated CP observable that reflects the changes in distance to the satellite. The CP measurement noise is between 0.3 and 3 millimeters [2].
2.2. Single Difference (SD)
For receivers at stations and , observing the same satellite s at times tr1 and tr2, one can write two PR equations like (2.1) and another two for the CP (2.3). Since and are approximately equal, the single difference (SD), illustrated in Figure 1, represents the difference between the observables collected from the two receivers and .
In the single difference equation, the errors related to the satellite (like orbit and satellite clock errors) are assumed the same for the observations collected from and . They are therefore cancelled out. Furthermore, if the baseline is short, in general, the ionospheric and tropospheric effects are spatially correlated and are also eliminated. The PR and CP SDs can thus be written as
where ∆ represents the difference between receivers, so , , , , and .
As the multipath depends on the geometry and the environment of each station, it is not eliminated in the SDs.
2.3. Double Difference
If two receivers and observe two satellites and at the same time, it is possible to form the double difference (DD) observable, that is, the difference between two SDs. The DD is illustrated in Figure 2.
The equations of PR and CP DDs can be written as
where ∇ represents the difference between the satellites, so , , , and .
An important feature of the DD is that the receiver clock errors ( and ) and the initial phases are eliminated. Therefore, when the DD is used in the functional mathematical model for short baselines and “similar’’ conditions at both stations, only the multipath effect is not eliminated. In the next section, we describe the proposed method to reduce the multipath from the DD and improve this functional model.
3. Multipath Mitigation Using Wavelets
The classical KF processing of GPS measurements generates residuals, which contain the signature of both nonmodeled errors and random measurement noise. If the functional mathematical model is adequate, the residuals obtained from the KF solution should be randomly distributed [11, 12]. However, as seen in the previous section, the multipath effect is not reduced in the functional model. It is desirable and important to extract the remaining errors from the residuals. Due to a lack of ways to do this, the multipath effect is generally not considered, mainly in kinematic positioning.
The detection and correction of the multipath from DD residuals using wavelet spectral analysis were described for static applications by [14, 15]. In those papers, the basis of the wavelet method was described. Here, the focus is mainly on presenting aspects necessary to apply the proposed wavelet method in kinematic positioning.
The multipath detection and mitigation using wavelet spectral analysis is based on four steps.
(1) Wavelet decomposition. Considering a relative kinematic positioning involving a short baseline, the DD residuals from the KF are decomposed by the wavelet transform using the pyramidal algorithm [10]. In this step, the Spline mother wavelet with five vanishing moments is used [15].
(2) Modify the wavelet coefficients by thresholding. The idea here is to detect and extract the multipath. In this sense, considering that the wavelet coefficients at the finest scale are, with a few rare exceptions, essentially pure noise, a nonlinear thresholding method [16] is applied at this level. The idea is to set to zero, by a hard threshold, those wavelet coefficients that do not exceed a suitable threshold λ. One approach commonly used for λ is the “universal” threshold . This is applied by the motivation that is asymptotically at least as smooth as the time series [17]. The standard deviation of the wavelet coefficients must be estimated. In this approach, suggested by [16], we used , with being the finest scale. Usually, for regression with nonstationary errors, a level-dependent λ becomes necessary [17]. But in this proposed method, nonstationary errors should not be excluded (set to zero). They need to be kept in the time series that represents the multipath errors. Therefore, just the undesired parts (noise) are removed and the remaining coefficients carry the significant information, which in this case (DD residuals from short baselines), is relative to the multipath components.
(3) Multipath correction. The multipath components from step 2 are then used for the reconstruction by the inverse wavelet transform (inverse pyramidal algorithm). Then, the extracted multipath errors are directly applied to the DD observables to correct them for the multipath effect.
(4) Improving the functional model. Once the DD observations are corrected for the multipath effect, the functional model (2.5) is improved:
At this step the KF is performed again, and the remaining residuals are now characterized just by the measurement noise.
Unlike in static applications, in kinematic multipath detection it is necessary to use data windows to process the data. Therefore, the four previous steps are performed for each window. Regarding step 1, in each instant k, the PR and CP DD residuals from instant to instant k are decomposed to estimate the multipath. As the pyramidal algorithm is used to apply the wavelet transform, it is necessary to consider . Thus, to process data in each instant, data from some previous instants are used in the multipath detection, but just the data from the current instant are corrected.
4. Analyses of the Functional Model in Kinematic Positioning
In this section, the DD functional model is analyzed in order to get insight on what it actually represents in practice. It is shown how the DD functional model, which involves differences of distances between receivers and satellites, describes the movement of the rover receiver. Therefore, only using DD observations time series one can have information about the rover movement.
To simplify the problem, we use graphic operations with vectors for a generic case. However, the process can be realized algebraically by finding the components of each vector, combining them to form the components of the resultant vector, and then converting to polar form.
The sum of two vectors, A and B, is a vector C, which is obtained by placing the initial point of B on the final point of A, and then drawing a line from the initial point of A to the final point of B. Vector subtraction is defined by the difference of two vectors generating a vector or . Thus, vector subtraction can be represented as a vector addition of the vector A with the vector B in the opposite direction. This process is illustrated in Figure 3 for an SD. Each vector represents a PR observable with a direction and magnitude, which is the distance between the receiver and the satellite.
(a)
(b)
(c)
Independent of the receiver and satellite positions, it is possible to take a plane passing by three points, similarly to that in Figure 3(c). For each SD formed, the resultant vector theoretically coincides with the baseline vector, that is, the distance between the receivers and . The several SD resultant vectors (C vector in Figure 3(c)) would be the same if the receiver clock and other errors were not present in the observables.
When a DD is computed, a difference is performed between two SD vectors, as illustrated in Figure 4. The small variations are shown in Figure 4.
We note that the proceeding refers to an analysis of the DD as a vector. However, the DD is used as a scalar to calculate the coordinates of the position. In this case, the differences of distances are realized, allowing that the geometry of the satellite distribution be taken into account; that is, the value of each SD depends on the position of each satellite.
In relative kinematic positioning, the receiver , called the base station, remains static while the rover is in movement. When the navigation starts, the DD vector direction depends on the movement direction of the rover receiver. Furthermore, the DD magnitude (vector length) is modified proportionally to the distance of the rover from the base station. If the rover is moving in a circle, the DD vector will also describe a circle trajectory. Figure 5 illustrates the DD vector in the instants ti and ti+k.
Some experiments, presented in the next section, confirm this theoretical analysis.
5. Experiment Description
An experiment was conducted at the “Luiz de Queiroz” College of Agriculture of São Paulo University, Piracicaba, Brazil, in November 9, 2007. This experiment was carried out in a kinematic mode with a controlled vehicle movement. The vehicle (small tractor) moved anchored around a pivot (Figure 6). The pivot has a bearing fixed in a reinforced concrete base. The tractor is fixed by a steel cable on a spin axis with minimum attrition. The receivers used, either in the pivot as a base station or in the rover station (tractor), were Topcon Hipper GGD.
In order to cause the multipath, a bus was placed near (2 m) the receiver (Figure 7). The bus remained stopped during the experiment.
The data were collected at a sample rate of 1 per second. The baseline reference length (14.98 m) was computed by topographic techniques, that is, by independent means (not GPS). This allows us to evaluate the accuracy of the proposed method in relation to the reference value.
Since the baseline length is short, errors resulting from ionosphere, troposphere, and orbits are assumed to be insignificant. Therefore, the DD estimated residuals for the PR and CP should exhibit mainly multipath and observable noise.
The L1 PR and CP DDs were processed using the software GPSeq, which is under development at UNESP, Presidente Prudente [18]. The wavelet method used to mitigate the multipath effect was also implemented in this software.
At the beginning of the experiment, the tractor remained stopped near the bus for about 20 minutes in order to solve the integer ambiguities. Then, the tractor started to move around the pivot. The static (pivot) antenna also rotates with respect to an axis through its center (spinning axis), thus experiencing the same phase wind-up effect as the slave. This common error is removed after performing an SD. The repeated loops provide ways of testing the results more than once.
PRN 14 with the highest elevation angle, as illustrated in Table 1, was chosen as base satellite to form the DD.
The results and analyses are presented in the next section.
6. Results and Analyses
The kinematic relative positioning was performed using two strategies. In the first one, hereafter referred to as the Standard, no multipath mitigation method was applied in the KF. In the second, the wavelet method (WAV) described in Section 3 was applied.
Not all collected data were processed in one run due to implementation restrictions. The processing is shown for a period of 1000 seconds of data, where the tractor remained stationary during the first 400 seconds and then started its movement, performing about 5 loops. Thus, it is possible to see the performance of the method for a static and kinematic processing as well as the transition between both cases. For other datasets, the results were similar.
In each instant k of the processing, the low-frequency multipath is detected, extracted, and directly applied to the DD observables to correct them for this effect. The KF is immediately performed again but with the corrected DD observables.
The first analysis performed is related to the aspects discussed in Section 4. Then, the improvement obtained by the WAV method is verified and presented in Section 6.2.
6.1. Functional Model Analysis
Figures 8, 9, 10, 11, and 12 compare the DD functional model behavior before and after applying the wavelet method for PR.
Figures 8 to 12 show the principles discussed in Section 4. One can observe the movement of the tractor, which started at about instant 400. When the navigation starts, the DD functional model represents the cyclic movement of the rover trajectory. It is also possible to observe that the WAV DD corrected for the multipath effects appears to better represent the movement of the rover receiver.
For the CP observables, the DDs are plotted in Figures 13, 14, 15, 16, and 17.
As the measurement noise for the CP is much smaller than for the PR, it is difficult to see differences between the Standard and WAV processing. But for either the PR or the CP DDs, it was possible to verify that the DD functional model represents the rover receiver movement.
6.2. Functional Model Improvement
To analyze the improvement in the functional model, the DD residuals obtained from Standard and WAV results were also compared. The estimated PR residuals are plotted in Figures 18, 19, 20, 21, and 22.
We note that the DD residuals involving PRN 18 with the lowest (55°–59°) elevation angles present the largest errors (Figure 20). On the other hand, the DD residuals involving PRN 03 with the highest elevation angles (23°–18°) present the smallest errors (Figure 19). These results agree with the fact that the multipath errors are inversely proportional to the satellite elevation angle.
In the same way, DD residuals obtained from the Standard and WAV strategies were also compared for CP residuals in Figures 23, 24, 25, 26, and 27.
One can observe from Figures 18 to 27 that the trend of the multipath effect was significantly mitigated for both the PR and the CP when the wavelet method was applied. However, the results for PR were much better, as is shown in the next analyses, where some statistics are presented. The Root Mean Square (RMS) for the DD residuals can be defined by [19]
where Vi represents the DD residual vector, and n is the number of observations used in the processing. The RMS and maximum error in the residuals for each DD are shown in Tables 2 and 3 for the PR and the CP, respectively. The reduction factor is computed by /.
Tables 2 and 3 verify that the improvement was better for the PR observable. This was expected because the multipath error is much larger for PR than for the CP. As a result, the correction for PR should be more significant.
To compare the quality of the PR and CP DD observables, with and without the multipath correction, the Local Overall Model (LOM) test statistic was applied. This test statistic has a Chi-squared distribution and allows a local detection of errors in each instant k. The definition of this test is given by [20]
where is the residual covariance matrix, and q is the degree of freedom in each instant k. The smaller the test value is, the better the quality of the observations is. The LOM test values are illustrated in Figure 28.
Figure 28 indicates that the LOM statistics improved significantly after using the wavelet method to correct the low-frequency multipath, showing that this effect was significantly mitigated. The LOM values from standard processing reached up to 38, while after the multipath mitigation, the maximum LOM value was about 4. The average improvement was about 80%.
In relation to the ambiguities solution, they could not be solved in the Standard processing because of the multipath effect. After correcting for this effect by the WAV method, the ambiguities could be reliably fixed at the beginning of the processing (less than 7 minutes).
In order to evaluate the accuracy of the baseline estimation, the baseline length estimated for each instant was compared with the reference value. For this, the discrepancies between the baseline estimative and the known baseline value were computed for standard and WAV processing. The RMS of the discrepancies for standard processing was 0.296 m. After applying the WAV method, the RMS was reduced to 0.029 m, representing a 90% improvement. Thus, one can verify that after the multipath correction, the results were better than in the standard processing.
7. Conclusion
In this paper, we show how to construct a functional mathematical model in order to eliminate most GPS errors.
The behavior of the functional model in kinematic mode is also demonstrated and analyzed in terms of vector operations. The experiment verifies the discussed principles, concluding that the DD functional model, which involves the differences in distances between receivers and satellites, describes the movement of the rover receiver.
We also present how the multipath effect, the only significant error not eliminated in the functional model, can be corrected by a proposed method using wavelets in kinematic positioning.
The results obtained in the experiments were very promising, showing that the proposed wavelet method appears to be a powerful method for mitigating multipath effects in kinematic GPS applications.
The multipath trend in the DD residuals was significantly corrected. The LOM statistical test also showed improvements in the quality of the data once the wavelet method was applied.
It was possible to observe the performance of the method for static and kinematic processing as well as the transition between both cases. The method demonstrated the potential for adaptation to the experimental condition. Other challenges must be addressed, such as more reflectors and different distances from the antenna, but due to the principles and features of the wavelet method, these can probably be accommodated.
Furthermore, the experiment was carried out with the rover varying smoothly. The next step is to evaluate how the solution is affected if the motion of the rover is highly dynamic. Free movements can also be experimented with, but the circular and controlled movement shown in this paper is very important for correct evaluation of the results.
The data processing and analyses have been restricted to the L1 observables but can also be extended to jointly process L1 and L2 observables.
Acknowledgments
The authors thank FAPESP for financial support to the first author (03/12770-3) and the Thematic Project (06/04008-2). They also thank Dr. José Paulo Molin and Thiago Martins Machado (“Luiz de Queiroz’’ College of Agriculture of São Paulo University) who made the kinematic experiment possible.