A Novel Method for Mitigating the GPS Multipath E ﬀ ect Based on a Multi-Point Hemispherical Grid Model

: The multipath e ﬀ ect is a crucial error source caused by the environment around the station and cannot be eliminated or mitigated by di ﬀ erential algorithms. Theoretically, the maximum value for the carrier phase is a quarter the wavelength, i.e., about 4.8 cm for the GPS L1 signal. Considering the increasing demands of high-precision applications, the multipath error has become a major factor a ﬀ ecting the accuracy and reliability of GPS millimeter-level data processing. This paper proposes a multi-point hemispherical grid model (MHGM) to mitigate the multipath e ﬀ ect. In this method, the hemisphere centered on each station is divided into a grid, and the multipath error at the station is estimated based on the parameterization of the grid points. The double-di ﬀ erenced (DD) observed-minus-calculated (OMC) values on some previous days are treated as the observation values to model the present multipath error. Contrary to the present methods which rely much on the platform of data collection and processing, MHGM can be potentially applied to GPS data processing with the existing hardware and software. Experiments in high-multipath and low-multipath environments are designed by mounting a ba ﬄ e or not. The experimental results show that MHGM is e ﬀ ective in mitigating the multipath e ﬀ ect. When using data from the previous day, an average improvement of about 63.3% in the RMS of DD OMC can be made compared with that without correction, and this is basically consistent with the sidereal ﬁltering (SF) method which is 63.0%. Furthermore, the e ﬀ ectiveness of the above two methods is better than that of the empirical site model (ESM). The kinematic positioning results are also basically consistent with the statistical results of the RMS values of DD OMC. Historical data from more than one day can more explicitly and e ﬀ ectively model the MHGM. Furthermore, compared with the SF method, the MHGM can be used not only to mitigate the multipath error, but also to orientate the sources of the multipath error around the station, and give guidance in the physical elimination of these sources.


Introduction
When the direct signals, reflected signals and diffracted signals of global positioning system (GPS) satellites reach the receiver antenna simultaneously, the interference will lead to an important error called the multipath error [1]. For carrier phase observations, the maximum error can be as residuals [32]. Ragheb (2007) compared SF methods based on the coordinate domain with those based on the observation domain, concluding that the former was more computationally efficient, but the latter was more accurate [33]. Zhong (2010), using SF on a short baseline of a single-differenced observation, was able to obtain more accurate positioning results than those possible with corrections applied in the coordinate domain [34].
As with modeling the multipath error in the coordinate or observation domain, based on the repetitive characteristics of satellite orbits, modeling the error in the spatial domain can also be an effective solution for mitigation of its influence. Cohen and Parkinson (1991) used observed data to set up a lookup table mapping the multipath environment around the receiver antenna [35]. Fuhrmann (2015) pointed out that it is unreasonable that stacking cells have a fixed azimuthal resolution, so an advanced method that makes use of congruent cells which have a similar shape and size was proposed [36]. Further, in its realization, the precise point positioning (PPP) residuals are used to generate the multipath stacking maps, which may add satellite clock errors and orbital errors in the final model. Moore (2014) used an empirical site model (ESM) to mitigate unmodeled site-specific errors [37]. Zero-differenced (ZD) carrier phase residuals are used to generate a site-specific correction map in this method. However, the technique to transform double-differenced residuals into ZD residuals may be not suitable due to the assumption of "zero mean" [38]. Dong (2016) proposed a multipath error correction method based on a multipath hemispherical map [39]. This method averages the residual of observed values corresponding to the fixed solution of single-difference ambiguity between stations of the satellite when a common-view satellite passes through a grid on the hemisphere between stations. Then, the low frequency multipath effect is retained and can be modeled while the noise effect is weakened. It solves the disadvantages due to the assumption of "zero mean", but its available distance is limited by a common receiver clock with multiple receivers.
In this work, we first give a brief presentation of the characteristics of the multipath effect in the space domain. Based on that, a novel method called the multi-point hemispherical grid model (MHGM) for mitigating the influence of the multipath error is introduced. The new method divides the hemisphere centered on each station into a grid, and estimates the multipath error at the station based on the parameterization of the grid points. Conventional double-differenced observations are used for modeling, so it can be potentially applied to GPS data processing with the existing hardware and software. Afterward, experiments in different multipath environments are designed, and observation residuals and positioning results before and after correction using MHGM are analyzed and compared with the results of the SF method and ESM. Finally, we draw some conclusions and give some suggestions of potential applications of MHGM in the future.

The Multi-Point Hemispherical Grid Model
In order to model the carrier phase multipath in the spatial domain, it is necessary to understand the characteristic of the multipath error. A model for the effects of the carrier phase multipath is built. Assuming that the receiver antenna is in a single planar horizontal reflector which is infinitely large and that its transmissions are only affected by a single reflected signal, the error M caused by the multipath effect on the carrier phase observation can be described by the following formula [40]: where H is the antenna height, ε is the angle of incidence of the reflected signal, α is the reflection coefficient of the reflected signal, and λ is the wavelength of the satellite signal. It can be seen from (1) that for a satellite signal of a certain wavelength, the multipath error depends upon the reflection coefficient, the angle of incidence of the reflected signal, and the spatial relationship between the antenna and the reflector. If the spatial relationship between the antenna and the surrounding environment remains unchanged during the observation process, changes in the Remote Sens. 2020, 12, 3045 4 of 18 propagation path of the reflected signal will mainly be due to changes in the direction of propagation of the satellite signal; its angle of incidence and reflection coefficient will change accordingly, ultimately affecting the multipath error present in the observed signal.
Based on the characteristics of the multipath effect in the space domain described above, a hemisphere can be defined centered upon the antenna's position, allowing the position of points to be described by a pair of angular coordinates (i.e., elevation and azimuth angles). Assuming an insignificant change in the environment around the station, and a given signal frequency, for a point in the hemispherical coordinate system based on the antenna, the multipath effect is related only to the elevation and azimuth angles of the satellite, and is not dependent on the signal's time of observation, or the satellite's pseudo-random noise (PRN) code.
This paper thus proposes the multi-point hemispherical grid model. This method aims to isolate multipath information, using inter-station double-differenced (DD) observed-minus-calculated (OMC), allowing the establishment of a multi-point hemispherical grid model consistent across stations. The lower hardware requirements of the novel approach allow it to be potentially applied in existing GPS network data processing systems. In particular, the use of DD OMC to model the multipath error does not require specialized receiver equipment, and the distance between stations is not limited by hardware considerations.
A grid is constructed for each station for the MHGM. For the elevation angle, the minimum value for grid points is set to E 0 , and the maximum value is set to E 1 . The azimuth angle of grid points is set to 0 • to 360 • , the range of possible values for it. A grid point is also set at the zenith of the coordinate system. Grid points are separated by an interval of d e for the elevation angle and d a for the azimuth angle. The number of parameters to be estimated in the model is determined by the value of d e and d a . The more parameters to be estimated, the more memory resources of the computer are consumed and the longer computation time is taken, but the description of the multipath error at the station is more detailed. The above factors should be taken into account to determine the value of d e and d a in modeling. An example of how the grid could be divided for a set of stations is given below (E 0 = 0 • , E 1 = 60 • , d e = 30 • , d a = 30 • ); the corresponding model parameters must be estimated for each of the black grid points shown.
The Positioning and Navigation Data Analyst (PANDA) software was used for attaining DD OMC between satellites and stations [41]. PANDA works under the ZD processing mode and the zero-differenced observed-minus-calculated (OMC) between each station and each visible satellite can be attained directly [42]. After that, the ZD ambiguities are mapped to a maximum set of independent DD ambiguities. Finally, the DD ambiguities could be solved and fixed sequentially [43]. For each set of fixed DD ambiguities processed by the PANDA software, the four groups of ZD OMC in the corresponding period could be obtained, so that the DD OMC of the carrier phase observation corresponding to the DD ambiguity is calculated according to the following equation as The ∆∇N in Equation (2) stands for the DD ambiguities. Normal equations corresponding to the grid points to be evaluated are established for the DD OMC s between the stations (m, n) and the satellites ( j, k), corresponding to an arbitrary set of fixed-ambiguity solutions. Figure 1 shows that the multipath error between the station m and the satellite j is assumed to involve four points Q  If the spatial relationship between the pierce point of the propagation of the signal of satellite and the hemispherical grid point parameters is as shown in the diagram on the upper right side of Figure 1, and it is assumed that the elevation angles and azimuth angles corresponding to the four grid point parameters 1 , 2 , 3 and 4 are respectively ( 1 , 1 ), ( 1 , 2 ), ( 2 , 2 ) and ( 2 , 1 ), then estimates of the value of model parameters at can be calculated by bilinear interpolation, as follows: .
If the spatial relationship between the pierce point of the propagation of the signal of satellite , and the hemispherical grid point parameters is as shown in the diagram on the lower right side of Figure 1, and assuming that corresponding elevation angles and azimuth angles of the three grid point parameters 1 , 2 and 3 are respectively ( 1 , 1 ), ( 2 , 2 ) and ( 2 , 3 ), then estimates of the value of model parameters at can be calculated by plane interpolation, as follows: .
Similarly, expressions for the values , and of model parameters, between station and satellite , station and satellite and station and satellite , can be established. Then, the DD OMC also can be expressed as follows: Observation (5) is used to derive normal equations corresponding to the relevant grid point parameters, and the fused solution, for several different pairs of satellites and ground stations and for data observed during different periods, is generated through superposition of these normal equations. The accuracy and reliability of the model's solution can also be further enhanced by imposing constraints on the correlation among grid point parameters, and imposing a requirement If the spatial relationship between the pierce point of the propagation of the signal of satellite j, and the hemispherical grid point parameters is as shown in the diagram on the lower right side of Figure 1, and assuming that corresponding elevation angles and azimuth angles of the three grid point parameters Q j m1 , Q j m2 and Q j m3 are respectively (e 1 , a 1 ), (e 2 , a 2 ) and (e 2 , a 3 ), then estimates of the value of model parameters at V j m can be calculated by plane interpolation, as follows: .
Similarly, expressions for the values V k m , V j n and V k n of model parameters, between station m and satellite k, station n and satellite j and station n and satellite k, can be established. Then, the DD OMC s also can be expressed as follows: Observation (5) is used to derive normal equations corresponding to the relevant grid point parameters, and the fused solution, for several different pairs of satellites and ground stations and for data observed during different periods, is generated through superposition of these normal Remote Sens. 2020, 12, 3045 6 of 18 equations. The accuracy and reliability of the model's solution can also be further enhanced by imposing constraints on the correlation among grid point parameters, and imposing a requirement that the absolute values of grid point parameters must be less than 1/4 the wavelength. The specific process of the model parameter estimation is displayed in Figure 2.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 19 that the absolute values of grid point parameters must be less than 1/4 the wavelength. The specific process of the model parameter estimation is displayed in Figure 2.

Experimental Design
Experimental verification was carried out on the top floor of the academic experiment building at Wuhan University during days 226-233 in 2018. Three GPS stations, A, B and C, were set up, with a metal baffle mounted to the west of station A to create a simulated high-multipath environment. Stations B and C were positioned in normal observation environments. The hardware specifications of stations A, B and C used in the experiment are shown in Table 1. The antenna used in each station was a common survey antenna without any anti-multipath techniques. The distances between the three stations were 2.7, 4.7 and 6.9 m, respectively, and their spatial relationship with the surrounding observation environments is shown in Figure 3.

Experimental Design
Experimental verification was carried out on the top floor of the academic experiment building at Wuhan University during days 226-233 in 2018. Three GPS stations, A, B and C, were set up, with a metal baffle mounted to the west of station A to create a simulated high-multipath environment. Stations B and C were positioned in normal observation environments. The hardware specifications of stations A, B and C used in the experiment are shown in Table 1. The antenna used in each station was a common survey antenna without any anti-multipath techniques. The distances between the three stations were 2.7, 4.7 and 6.9 m, respectively, and their spatial relationship with the surrounding observation environments is shown in Figure 3.  Using the accurately known coordinates of the three stations, multipath error modeling was applied to data from L1 observations on days 226-233 in 2018. The effectiveness of multipath error correction for the observations from the 233rd day of 2018 was analyzed. The models and parameters used in the data processing of the PANDA software are listed in Table 2. The maximum DD multipath error should not exceed ½ the wavelength when the multipath errors reach the theoretical maximum and the signs are opposite at both ends of the baseline. For this reason, the DD OMC larger than 1/2 the wavelength will be eliminated in the data preprocessing before the process of modeling.  Using the accurately known coordinates of the three stations, multipath error modeling was applied to data from L1 observations on days 226-233 in 2018. The effectiveness of multipath error correction for the observations from the 233rd day of 2018 was analyzed. The models and parameters used in the data processing of the PANDA software are listed in Table 2. The maximum DD multipath error should not exceed 1/2 the wavelength when the multipath errors reach the theoretical maximum and the signs are opposite at both ends of the baseline. For this reason, the DD OMC larger than 1/2 the wavelength will be eliminated in the data preprocessing before the process of modeling.  For each day, the same DD ambiguities fixed results are used to establish MHGM, ESM and SF models. Among them, the MHGM also has been modeled using multi-day historical data. Therefore, this paper provides the results with different kinds of multipath error correction and the result without correction for comparison. The specific solution strategies adopted in the experiment (i.e., combinations of methods used and datasets to which they were applied) are shown in Table 3.

Analysis of Results
PANDA works under the ZD processing mode and the zero-differenced observed-minus-calculated between each station and each visible satellite can be attained directly. Corresponding to the strategies shown in Table 3, Figure 4 shows the ZD OMC for each satellite at the three stations A, B and C during a period of fixed ambiguity, in the form of a sub-satellite point track. With no correction for the multipath error (N), the results for all three stations, A (N), B (N) and C (N), include relatively large ZD OMCs. This is particularly noticeable for paths via the west side of station A where the baffle was mounted. So, the multipath effect is particularly obvious on the left side of diagram A (N) with many groups of large observation residuals. For the SF method of multipath error correction, the satellites' orbit repeat period advance must be taken into account. In this experiment, the SF in the observation domain is adopted, in which the average orbit repeat period of the whole GPS constellation is set to 23 h 55 m 55 s (86,155 s), with a daily advance of 245 s [28]. The results of A (SF), B (SF) and C (SF) show that the multipath effect has been effectively mitigated most of the time, but the residuals towards the end of the satellite trajectories are larger and it is a sequence that cannot be corrected effectively using the SF model. Using the MHGM proposed in this paper and the ESM, even only using the data of the previous day for modeling, as shown in A (M1/E1), B (M1/E1) and C (M1/E1), these uncorrected satellite trajectories are less obvious (i.e., the phenomenon of larger residuals towards the end of the satellite observation periods does not exist). The modeling effectiveness of MHGM can be further improved by using data observed during the prior two or more days. This is because the orbit repeat period of GPS satellites is about one day, and the multipath error in GPS observation data can be effectively modeled using data observed over the most recent two days. In order to verify the validity of MHGMs constructed by using multi-day GPS data, the M7 strategy is also shown in the following conclusions. Figure 4, abbreviated for reasons of space, shows only the ZD OMC at the three stations for the strategies of N, SF, ESM, M1 and M7.
To quantify the effects of the MHGM proposed in this paper in comparison with those of the SF and ESM methods, and to evaluate differences in the effectiveness of multipath error modeling based on data observed over 1~7 days, Table 4 shows the mean root mean square (RMS) of DD OMC during fixed-ambiguity periods on day 233 for all 10 strategies.
The HMP in Table 4 stands for the statistical results of baseline A-C which is in a high-multipath environment. If data observed on the previous day only are used for multipath error modeling, the effect of the ESM method is very obvious. The DD RMS improvement of strategy ESM compared with strategy N can reach 68.3%, but it is worse than that of SF. For strategy M1, it is comparable to the improvement level of strategy SF, and the improvement increases from 72.6% to 75.4%. For multipath error modeling based on data observed over the previous two days or more, the improvements offered by M2-7 average about 75.6%. The improvement between M2 and M7 is not very obvious, that is because MHGM can make full use of the correlation of the multipath error in the neighboring directions, and then the direction without valid modeling data also can be appropriately corrected by extrapolation or interpolation.   The LMP in Table 4 stands for the statistical results of baseline B-C which is in a low-multipath environment. Strategy M1 is slightly worse than that of strategy SF, but the results of multiple days can further improve the mitigation effect of MHGM. Figure 5 provides a graphical display of the distribution of the RMS values of the observation residuals from the double-difference ambiguity fixed arcs from the 233rd day, derived using the different processing strategies. As can be seen in Figure 6, the low-frequency patterns are basically eliminated after the correction using different models in the double-differenced residual series of G01 and G03 on station A and station B. Compared with SF, the methods based on the spatial domain can also mitigate the low-frequency part of the multipath effects, and the series are slightly more stable with model correction using MHGM. Power spectrum analysis indicates that the low-frequency multipath effects Distribution of RMS DD OMC for a fixed-ambiguity period under different processing strategies.
As can be seen in Figure 6, the low-frequency patterns are basically eliminated after the correction using different models in the double-differenced residual series of G01 and G03 on station A and station B. Compared with SF, the methods based on the spatial domain can also mitigate the low-frequency part of the multipath effects, and the series are slightly more stable with model correction using MHGM. Power spectrum analysis indicates that the low-frequency multipath effects can be reduced with these multipath corrections. In the frequency band higher than 0.1 Hz, however, the power of the MHGM is similar to the ESM and SF.
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 19 Figure 6. The double-differenced residual series and power spectral density in the period of fixed ambiguity, 233th day, 2018.
Precise positioning is a more intuitive form to assess the effect of multipath error correction. For this reason, the influence of multipath error correction on the positioning results using SF, ESM and MHGM methods is processed and shown in Table 5. Station C is fixed in the data processing, and the three-dimensional coordinates of static stations A and B are solved according to the kinematic relative positioning mode. For the new method, the strategies M1 and M7 are used to model the multipath error.

RMS (mm) A (HMP) B (LMP) Horizontal
Vertical Horizontal Vertical Precise positioning is a more intuitive form to assess the effect of multipath error correction. For this reason, the influence of multipath error correction on the positioning results using SF, ESM and MHGM methods is processed and shown in Table 5. Station C is fixed in the data processing, and the three-dimensional coordinates of static stations A and B are solved according to the kinematic relative positioning mode. For the new method, the strategies M1 and M7 are used to model the multipath error. When the multipath error is uncorrected, the positioning accuracy of station A is worse than that of station B because a metal baffle is mounted to station A to simulate the high-multipath environment.  Figure 7.
Compared with the SF method, the MHGM proposed in this paper not only enhances the improvements obtained through multipath error correction, but also effectively provides a graphical display of multipath error interference around the stations. Figure 8 shows the hemispherical modeling results from strategy M1 at the stations. In this, "overlook", "west" and "east" are visualizations of the effect from above, looking to the west side and looking to the east side, respectively. The top view shows higher modeled values in large areas of the left side area of A (M1). Further, in the west side view, an obvious and continuous abnormal area at the bottom is visible, the location of which coincides with the position where the metal baffles were mounted: this reflects the influence of these baffles on the multipath error at station A. This characteristic of the MHGM implies it can not only be used to mitigate the influence of the multipath error at stations, but also to orientate the sources of the multipath error around the station, providing guidance as to the measures which could be useful for physically eliminating the influence of these error sources.
In order to show the changes in the model results over multiple days intuitively, the modeling results of the three stations based on daily GPS observations from 1 day to 7 days are also attached to the Supplementary Materials (Supplemental-Gif-1.mp4). Due to space considerations, only the modeling results of MHGM using seven daily observations are shown below. It can be seen that the modeling results of MHGM from strategy M7 are generally similar to those of strategy M1, but the results of strategy M7 have more details in the top view of Figure 9. As previously mentioned, the orbit repeat period of GPS satellites is about one day, so multipath error modeling based on one-day observations will be feasible, and modeling base on multi-day observations can further make up for the lack of data that may exist in one-day observations. kinematic positioning results are basically consistent with the statistical results of the RMS values of DD OMC. In order to further demonstrate the influence of five processing strategies, N, SF, ESM, M1 and M7, on the positioning results at different times, the real-time kinematic positioning error sequences and the corresponding RMSs of station A (HMP) in the north, east and up directions are given in Figure 7.  Further, in the west side view, an obvious and continuous abnormal area at the bottom is visible, the location of which coincides with the position where the metal baffles were mounted: this reflects the influence of these baffles on the multipath error at station A. This characteristic of the MHGM implies it can not only be used to mitigate the influence of the multipath error at stations, but also to orientate the sources of the multipath error around the station, providing guidance as to the measures which could be useful for physically eliminating the influence of these error sources. In order to show the changes in the model results over multiple days intuitively, the modeling results of the three stations based on daily GPS observations from 1 day to 7 days are also attached to the Supplementary Material (Supplemental-Gif-1.mp4). Due to space considerations, only the modeling results of MHGM using seven daily observations are shown below. It can be seen that the modeling results of MHGM from strategy M7 are generally similar to those of strategy M1, but the results of strategy M7 have more details in the top view of Figure 9. As previously mentioned, the orbit repeat period of GPS satellites is about one day, so multipath error modeling based on one-day observations will be feasible, and modeling base on multi-day observations can further make up for the lack of data that may exist in one-day observations.

Conclusions
This paper proposes a novel method for the mitigation of multipath error at stations based on a multi-point hemispherical grid model. The method divides the hemisphere centered on each station into meshes, and estimates the multipath error at the station based on the parameterization of the grid points. During the process of grid point parameter estimation, the multi-point hemispherical

Conclusions
This paper proposes a novel method for the mitigation of multipath error at stations based on a multi-point hemispherical grid model. The method divides the hemisphere centered on each station into meshes, and estimates the multipath error at the station based on the parameterization of the grid points. During the process of grid point parameter estimation, the multi-point hemispherical grid model for each station is obtained, taking the DD OMC corresponding to the solutions of fixed ambiguity for different satellite pairs and for different periods as input information.
The experimental results show that for the multipath error modeled by data from only the previous day, the correction performance of the MHGM can be generally consistent with the SF method, and better than ESM. When the previous two or more days' observation data are used for modeling the multipath error, the MHGM can be described in more detail. The average improvement of the mean RMS of DD OMC increases from 63.3% (M1) to 66.3% (M7), in comparison to the uncorrected results. Further, the kinematic positioning results also coincide with the statistical RMS values of DD OMC. The results also show that the multipath error modeling from one-day observations will be feasible and modeling from multi-day observations can further make up the lack of data that may exist in one-day observations.
Compared with the SF method, the proposed MHGM can not only mitigate the multipath error effectively, but also provide a graphical distribution of multipath error interference around the stations, and thus give guidance in physically eliminating the sources of the multipath error in the future.
Author Contributions: Y.W., X.Z. and W.T. provided the initial idea and designed the research; X.Z. and Y.W. processed data, and wrote the manuscript; Y.L. helped to accomplish some tests; C.D., Y.Z. and J.F. helped with the writing. All authors have read and agreed to the published version of the manuscript.