Multipath mitigation for GPS/ Galileo/BDS-3 precise point positioning with overlap-frequency signals

The multipath effect is a major Global Navigation Satellite System (GNSS) error source due to its environment-dependent characteristic, which complicates its mitigation process for the high-rate determination of displacements. For instance, Sidereal Filtering (SF) and Multipath Hemispherical Map (MHM) require the observations spanning at least one full cycle of satellite orbit repeat period ( e.g. , ten days for Galileo navigation satellite system (Galileo) to reproduce the satellite geometry against ground stations. As a consequence, the practicability of SF and MHM is limited due to potential station-surrounding changes over a long period. In this study, we used the overlap-fre-quency signals on Global Positioning System (GPS) L1/L5, Galileo E1/E5a, and BeiDou-3 Navigation Satellite System (BDS-3) B1C/B2a to construct an interoperable MHM ( i.e. , MHM_GEC) across constellations to mitigate multipath more efficiently. We thus used 31 days of 1-Hz GPS/Galileo/BDS-3 data at 21 stations in Europe to compare this overlap-frequency MHM with those GNSS-specific MHMs ( i.e. , MHM_G for GPS, MHM_E for Galileo, and MHM_C for BDS-3), as well as SF. It is confirmed that the multipath effects on overlap-frequency signals are of a high spatial consistency across all GNSS. The mean reduction rate of applying MHM_GEC to GPS, Galileo, and BDS-3 carrier-phase residuals is 25%, 31%, and 28.5%, respectively, which are up to 25 percentage points higher than those of MHM_G, MHM_E, and MHM _C. Furthermore, the MHM_GEC constructed using 5 to 6 days of data can improve the positioning precision by 40%, outperforming the MHM_E, MHM_C, and SF using 10 days of data. Therefore, the interoperable MHM_ GEC is more efficient in mitigating multipath effects for high-precision GNSS positioning.


Introduction
High-rate Global Navigation Satellite System (GNSS) positioning has gained comprehensive applications in seismology, seismic engineering, meteorology, and other fields (e.g., Gang et al., 2013;Kouba, 2003;Larson et al., 2003).Unlike traditional inertial seismic methods, highrate GNSS positioning can directly capture instantaneous ground displacements, including low-frequency motions and even static offset (Blewitt et al., 2016;Bock et al., 2011).For instance, during the caldera collapse at Kīlauea Volcano in 2018, high-rate GNSS observations recorded the ground deformation in unprecedented detail (Anderson & Johanson, 2022).Unfortunately, high-rate GNSS determined displacements can usually only achieve centimeter-level precision, which will mask typical shortterm deformation signals of interest (e.g., Wright et al., 2012).The major source of error is multipath effects, which are attributed to satellite signal reflections by the surroundings of ground antennas (Wang et al., 2015).These non-line-of-sight signals can result in pseudorange multipath over tens of meters and carrier-phase multipath up to a quarter of a cycle (Dong et al., 2016;Moradi et al., 2015;Ragheb et al., 2007).In addition, GNSS multipath effect depends on numerous stationrelated factors (Elósegui et al., 1995) and signal frequencies (Georgiadou & Kleusberg, 1988;Lau & Cross, 2007).As a result, multipath effect is difficult to model accurately using mathematical equations or eliminate using observation differencing techniques (Lau & Cross, 2007;Sun et al., 2019;Xu et al., 2020).
At (quasi-)static GNSS stations, the method of correcting target GNSS observations or displacements based on multipath temporal and spatial repeatability has been proven effective.Genrich and Bock (1992) proposed Sidereal Filtering (SF), which capitalizes on the fact that the geometry of the Global Positioning System (GPS) satellite-reflector almost repeats every sidereal day, allowing the data from the preceding sidereal day to be used to correct the target observations.However, the repeated orbital period of satellites may deviate from theoretical predictions by several seconds to tens of seconds (Choi et al., 2004;Shen et al., 2020;Su et al., 2018).To address this issue, Zhong et al. (2010) and Dong et al. (2016) developed Advanced Sidereal Filtering (ASF), which calculates the orbit repeat period of each satellite using broadcast ephemeris.However, for the GNSS constellations with a long orbit repeat period, such as BeiDou-3 Navigation Satellite System (BDS-3) Medium Earth orbits (MEOs) and Galileo navigation satellite system (Galileo) satellites, ASF may require the observations of up to 7 and 10 sidereal days, respectively (Geng et al., 2017(Geng et al., , 2018)).During the long observation period, station surrounding changes may compromise the effectiveness of ASF (Ye et al., 2015).In addition, ASF corrections are specific to individual satellites and epochs, suggesting that the correction calculated for a satellite can only be used for that specific satellite (Wang et al., 2018).
Unlike the ASF, the Multipath Hemisphere Map (MHM) relies on the spatial location of a satellite, such as the satellite's azimuth and elevation (Cohen & Parkinson, 1991;Dong et al., 2016).The method assumes that multipath effects only relate to the satellite location in space if the station surrounding does not change.Regular grids or spherical harmonics can therefore be used to model and quantify this effect in the hemispherical skyplot, with more detailed multipath characteristics if satellite tracks are densely and evenly distributed.To achieve this, more days of data can be used.Zheng et al. (2019) used 15 days of data from Bei-Dou-2 Navigation Satellite System (BDS-2) MEOs and 20 days of data from Galileo satellites.However, like ASF, such modeling requires long observation periods, and potential station surrounding changes over long periods may undermine the model validity.In addition, MHM is primarily used for static monitoring stations, but can also be applied to mobile platforms when the primary multipath effect stem from the platform body itself (Cai et al., 2016).
Another way to enhance the density of satellite tracks in the hemispherical skyplot is to involve more satellites in the MHM modeling.Currently, there are 32 GPS satellites, 26 Global'naya Navigatsionnaya Sputnikovaya Sistema (GLONASS) satellites, 28 Galileo satellites, and 50 BDS (including BDS-2 and BDS-3) satellites available.However, simply gathering all these satellites together does not make sense, since the observations from different constellations may have different multipath effects due to different signal frequencies.Therefore, most of research publications have the limitation in constructing GNSS-specific MHM (Lu et al., 2021;Ren et al., 2023;Tang et al., 2021;Wang et al., 2019;Yuan et al., 2023).Fortunately, the GPS L1/L5, Galileo E1/E5a, and BDS-3 B1C/B2a signals (hereinafter referred to as overlap-frequency signals) share the same carrier frequency, making it possible to integrate these signals to establish an interoperable MHM model (Lu et al., 2023), which means that BDS-3 multipath can be used to correct Galileo multipath, and vice versa.By combining GPS, Galileo, and BDS-3 satellites, we can have much more observations within a few days, resulting in more evenly and densely distributed satellite tracks in the hemispherical skyplot.This approach can even achieve a spatial resolution comparable to GNSSspecific MHMs constructed using 10 days of data.Furthermore, this strategy can also reduce the risk of degrading MHMs over long span of data, as is the case of ASF and GNSS-specific MHM.
It is promising to use overlap-frequency signals for interoperable MHM.However, several questions need to be answered.First, can these signals provide a positioning performance equivalent to legacy-frequency signals such as GPS L1/L2, Galileo E1/E5a, and BDS-3 B1I/B3I?Second, is the multipath effect of overlapfrequency signals consistent across GNSS?Finally, how effective is the overlap-frequency MHM in mitigating multipath effects?We will first compare the positioning precision of GPS, Galileo, and BDS-3 Precise Point Positioning (PPP) using legacy-frequency signals and overlap-frequency signals.Then, we will determine the consistency of multipath effects in overlap-frequency signals through experimental verification.Finally, we will establish an interoperable MHM using overlap-frequency signals and compare its correction performance with GNSS-specific MHMs and the ASF methods.

Precise point positioning
PPP is performed to obtain the carrier-phase residuals.We use the ionosphere-free combination to eliminate the impact of low order ionospheric delay, while high order ionospheric delay is generally negligible.The equation for ionosphere-free carrier-phase observation is as follows: where s is the satellite index, � (s)  IF is the ionosphere-free carrier-phase observation, and ρ (s) denotes the geomet- ric distance between the receiver and the satellite; c , dt r , dt (s) , IF , and N (s)  IF represent the speed of light, receiver clock errors, satellite clock errors, wavelength, and integer ambiguity of ionosphere-free combination, respectively; IF , and ε (s) IF denote troposphere effect, receiver hardware bias, satellite hardware bias, multipath effect, and the observation noise of ionosphere-free combination, respectively.
When the exact receiver positions are known, ρ (s) can be easily calculated.The receiver clock offset dt r is esti- mated assuming white noise and separately estimated for each GNSS to prevent the influence of Inter-System Biases (ISB).Since the GPS satellite clock products provided by the International GNSS Service (IGS) analysis center were computed using the L1/L2 ionosphere-free combination, this resulted in inter-frequency clock bias when using the L1/L5 ionosphere-free combination for PPP.Therefore, we estimated the clock products using the legacy-and overlap-frequency signals separately to correct the corresponding satellite clock errors dt (s) .The ambiguity N (s)  IF can be estimated as a float constant.Noted that no integer ambiguity resolution was performed in this study.Satellite and receiver hardware delays, i.e.B (s)  IF and B r,IF will be assimilated into ambigu- ity and clock parameters.For the tropospheric delay T (s)  IF , its wet and dry delays are corrected by the Saastamoinen model (Saastamoinen, 1972) and its residual zenith delays are estimated as piece-wise constants every hour using the Global Mapping Function (GMF).By subtracting the calculated terms from the left-hand side of Eq. 1, we can obtain the residual term where R (s)  IF is the carrier-phase residual of ionosphere-free combination, M (s)  IF and ε (s) IF denote the multipath effect and observation noise of ionosphere-free combination.

Constructing and correcting for MHMs
In this study, the carrier-phase multipath errors from GPS L1/L5, Galileo E1/E5a, and BDS-3 B1C/B2a ionospherefree combinations are utilized to establish the MHM.As a pre-processing, we employ a low-pass filter to truncate the high-frequency measurement noise of the carrier-phase residuals: where L() represents the low-pass filter with an empirical cut-off frequency of 0.1 Hz, and M s IF denote the filtered carrier-phase residual, which is presumed as multipath errors.
In this study, a grid size of 1° × 1° is chosen for the MHM.The multipath errors within a grid are averaged as the multipath value in the MHM model (Dong et al., 2016): where a z and e l are the initial azimuth and elevation of each grid; a (i)  z and e (i) l are the azimuth and elevation of the i-th multipath error within this grid; n represents the number of multipath errors within this grid.To safeguard the accuracy of multipath correction from the interference of outliers, any multipath error values within a grid that surpass 2 times the sigma will be identified as outliers and consequently excluded.The number of remaining multipath errors is required to be equal to or greater than 30, and a threshold is chosen to strike a balance between the precision and resolution of multipath correction.Otherwise, the multipath corrections will be set to zero.
The multipath correction can be obtained from the MHM based on the elevation and azimuth of the satellite relative to the station: where a s z and e s l are the azimuth and elevation of the satellite, which fills in the MHM grid ( [a z , a z + 1); [e l , e l + 1) ); C MHM (a (s)  z , e (s) l ) denotes the obtained multipath correc- tion that is subtracted directly from the carrier-phase ionosphere-free combination to compensate for the multipath delay in the kinematic PPP process.If there is no correction value for a specific grid, the multipath correction is zero.

Constructing and correcting for ASF
The ASF is specific exclusively to satellites and epochs and can be modeled as: where the multipath error M (s) IF (t) of the satellite s at epoch t is derived from Eq. 3. We use the daily broad- cast ephemeris to calculate the repetition period in each epoch.The formula for calculating the satellite orbit repetition period using Kepler's third law is: where t s is the orbit repetition period of the satellite, k is the number of satellite revolutions,

√
GM is the standard gravity parameter, a is the semi-major axis of the satel- lite orbit, and n is the correction for the satellite angular velocity.Due to the presence of multiple sets of parameters in a daily broadcast ephemeris, we average the calculated values for each set of parameters as the satellite orbital repetition period: where t (i)  s is the satellite orbit repetition period calculated by the i-th parameter group of the satellite in the broadcast ephemeris, and t m is the averaged orbit repetition period of the satellite for that day.
The multipath correction can be obtained by substituting the satellite s (c) and its corresponding epoch t (c) into the ASF model: Then, we subtract C ASF (s (c) , t (c) ) from the carrier-phase ionosphere-free combination to compensate for the multipath delay at each epoch in the kinematic PPP process.
Figure 1 illustrates the flow of multipath modeling and correction for overlap-frequency signals using MHM and ASF.We first utilize static PPP to calculate the residuals of the ionosphere-free combination of overlap-frequency signals from GPS/Galileo/BDS-3 (Eqs. 1 and 2).Then we apply low-pass filtering to extract the multipath delay from the residuals (Eq.3).The obtained multipath delays in conjunction with the satellite azimuth and elevation are used to create the MHM (Eq.4), and with the corresponding time to create the ASF model (Eq.6).By substituting the target day's satellite azimuth and elevation angles into the constructed MHM (Eq.5) and the target day's satellite orbital repetition period (Eqs.7 and 8) into the constructed ASF model (Eq.9), we can derive the corresponding kinematic PPP multipath corrections.( 6)

Data processing
For algorithm validation, 21 IGS stations in Europe were chosen.Among them, 9 peripheral ones marked in green were specifically utilized to estimate the 1 Hz satellite clock products, reducing the impact of orbital errors.The remaining 12 stations in various observation environments highlighted in red were used for PPP, as shown in Fig. 2. Further, the stations KOS1 and SAS2 marked with red pentagon stars have relatively low multipath effects due to their open observation environments.Conversely, the other ten stations marked with red circular dots suffer from severer multipath effects.All 21 stations can track the multi-frequency signals of GPS/Galileo/BDS-3 constellations, including overlap-frequency signals like L1/E1/B1C and L5/E5a/B2a, as well as legacy-frequency signals such as L1/E1/B1I and L2/ E5a/B3I.We utilized their 1-Hz observations of a total of 31 days from day 290 to day 320 of 2021.On average, there are about 4.7 GPS satellites available for overlapfrequency signals and 7.3 GPS satellites for legacy-frequency signals.Additionally, there are about 6.9 BDS-3 and 7.8 Galileo satellites usable at each epoch.The mean value of Position Dilution Of Precision (PDOP) is about 1.9 for overlap-frequency signals and 1.7 for legacy-frequency signals at each station.Table 1 lists each station's antenna, radome, receiver type, firmware version, and supported signal frequencies.
The last day of each segment serves as the target day for the multipath correction, i.e., days 300, 310, and 320.The observations from the target day backward to the first day of the segment, ranging from 1 to 10 days, were utilized to establish the MHM with different durations.For example, the 1-day MHM for day 300 is based on the data from day 299, the 2-day MHM is based on the data from days 298 and 299, and so on until the 10-day MHM.
We conducted 1-Hz clock estimation for legacy-frequency and overlap-frequency observations using the nine peripheral stations (green triangle in Fig. 2).Using such a regional network to estimate clock products can help absorb the common errors of the regional stations, particularly in reducing the impact of orbital errors.During the estimation process, we utilized the one-second interval measurements from these nine stations and applied the clock estimation method proposed by Geng et al., (2019aGeng et al., ( , 2019b)).The raw clock and orbit products were obtained from the European Space Agency (ESA), with a sampling rate of 30 s.With these 1-Hz clock products, we then performed static PPP with combined GPS/ Galileo/BDS-3 to derive float carrier-phase residuals.These carrier-phase residuals were then used to construct the MHM and ASF models, as illustrated in Fig. 1.Note that the station coordinates were fixed for both the clock products estimation and the static PPP.Ultimately,  (Geng et al., 2019a).We have developed a multipath modeling module for PRIDE PPP-AR version 3.1 based on the algorithm in this article (https:// github.com/ Pride Lab/ PRIDE-PPPAR).

High-rate positioning using GPS/Galileo/BDS-3 overlap-frequency signals
To assess the efficacy of positioning with overlap-frequency signals, specifically, to determine if it's feasible to achieve a positioning precision comparable to or surpassing that with legacy-frequency signals, we conducted a comparative analysis of positioning precision of PPP utilizing these two types of frequency signals individually.Figure 3 shows the positioning differences in the east, north, and up components of PPP when using overlap-and legacy-frequency signals, i.e., GPS L1/L5 and L1/L2, Galileo E1/E5a and E1/E5a, BDS-3 B1C/B2a     and B1I/B3I.Although the overlap-and legacy-frequency of Galileo are the same, there are slight differences in their PPP results due to the use of different re-estimated clock products.Due to the constrained availability of GPS satellites in the L5 frequency, we only calculated the results during the periods when the number of satellites exceeded 6 during GPS positioning.
The PPP differences using the overlap-frequency and legacy-frequency signals for all 12 stations over 31 days, as shown in Fig. 3, exhibit small variations.For GPS, the mean positioning differences in the east, north, and up components are merely 0.05 cm, 0.03 cm, and 0.19 cm, respectively.Furthermore, the percentages of positioning differences in these three components smaller than 0.5 cm, 0.5 cm, and 1 cm are 90%, 79%, and 79%, respectively.For Galileo, the positioning differences are almost negligible.This is because the two signals use the same frequency, and the differences are primarily due to the re-estimated clock products used.For BDS-3, the mean positioning differences in the east, north, and up components are 0.09 cm, 0.2 cm, and 0.19 cm, respectively, slightly larger than the other systems.However, the percentages of positioning differences in these three components smaller than 0.5 cm, 0.5 cm, and 1 cm remains high, 90%, 83%, and 80%, respectively, indicating that the differences are still small.These results demonstrate that the PPP with overlap-frequency signals can achieve the precision comparable to that of using legacy-frequency signals.
Table 4 compares the mean Root Mean Square (RMS) of PPP positioning errors for 12 stations over 31 days using the BDS-3 legacy-and overlap-frequency signals.
The results show that the positioning precision with the overlap-frequency signals is superior to that with the legacy-frequency signals at all 12 stations.Specifically, the mean RMS of positioning errors in the east, north, and up components with overlap-frequency signals at the 12 stations are 1.08 cm, 0.99 cm, and 2.06 cm, respectively.This represents a reduction of 15.6%, 18.3%, and 17.3% compared to that with the legacy-frequency signals which have RMS values of 1.20 cm, 1.21 cm, and 2.50 cm, respectively.The improved performance with overlapfrequency signals can be attributed to the lower pseudorange noise and higher received power of the B1C and B2a signals compared to the B1I and B3I signals (Geng et al., 2023).
Furthermore, Table 5 provides a comparison of the mean RMS of positioning errors using GPS legacy-and overlap-frequency signals at the same 12 stations over 31 days.The results show that the positioning precision of GPS overlap-frequency signals is slightly higher than that of legacy-frequency signals, but the improvement is not as significant as that of the BDS-3.Specifically, the mean RMS of positioning errors in the east, north, and up components with overlap-frequency signals at the 12 stations are 0.88 cm, 0.84 cm, and 1.87 cm, respectively.This represents a reduction of only 5%, 3%, and 4% compared to the legacy-frequency signals.This small improvement may be attributed to the larger bandwidth of the L5 signal used in the overlap-frequency signals, which provides lower observation noise compared to the L2 signal used in the legacy-frequency signals (Bolla & Borre, 2019).Overall, Tables 4 and 5 demonstrate the comparable or better positioning precision using overlapfrequency signals compared to legacy-frequency signals can be achieved.

Multipath consistency among GPS, Galileo, and BDS-3 on overlap-frequency signals
Figure 4 presents an illustrative example showcasing the color-coded carrier-phase multipath delays of the overlap-frequency signals of GPS, Galileo, and BDS-3 satellites as they traverse a particular 1° × 1° skyplot grid.The data spans 10 days and includes a total of 773 GPS epochs, 1008 Galileo epochs, and 1339 BDS-3 epochs.The mean multipath errors for all GPS, Galileo, and BDS-3 satellites are -0.018cycles, -0.019 cycles, and -0.018 cycles, respectively, with STDs of 0.011 cycles, 0.012 cycles, and 0.010 cycles.It can be observed that the orbits are mostly represented by the light blue color in this grid, with a STD of 0.011 cycles.Consequently, Fig. 4 provides a visual evidence supporting the notion that the GPS, Galileo, and BDS-3 satellites exhibit a consistency in terms of their overlap-frequency multipath effects within a specific skyplot grid.
Figure 5 further illustrates the STD of multipath errors for the overlap-frequency signals at the station BAUT in each 1° × 1° grid during the period from day 310 to day 319.To ensure statistical reliability, each grid contains at least 30 epochs (see Eq. 4).The mean STD for multipath errors across all grids is 0.0069 cycles.In Fig. 5, most of the grids are represented in blue, indicating that about 90% of the grids have STD of less than 0.01 cycles.Only grids with the elevation angles less than 15° have STD of more than 0.02 cycles, which are represented in purple.This is primarily due to larger variations in multipath errors within grids at low elevations.Furthermore, the situation for the other two 10-day segments (i.e., days 290-300 and days 300-310) is similar.Their mean STD of multipath errors across all grids is 0.0065 cycles and 0.0072 cycles, respectively.Almost all grids in these segments have STDs of less than 0.03 cycles.Therefore, Figs. 4 and 5 demonstrate the general spatial consistency of the overlap-frequency multipath effects of GPS, Galileo, and BDS-3.This consistency serves as a key basis for establishing an interoperable MHM.

Comparison of four MHM models' coverages
Figure 6 depicts the spatial distribution of one-day carrier-phase multipath corrections for overlapfrequency signals of the GPS, Galileo, and BDS-3 constellations.The color distributions exhibit striking similarities, such as the presence of a blue band in the central position (i.e.elevation 30-60° and azimuth 90-270°) across all panels.This indicates a consistent spatial distribution of multipath errors in overlap-frequency signals from different constellations.It is worth noting that significant blanks can be observed in the skyplot of the single-constellation (i.e.first three panels).By combining Galileo and BDS-3 or combining GPS, Galileo, and BDS-3, the blanks in the skyplot are noticeably reduced.Additionally, the similarity in the spatial distributions of the corresponding multipath errors becomes more evident.Consequently, by combining constellations, it is possible to leverage the same multipath spatial distribution characteristics of the overlap-frequency signals, thereby enhancing the completeness of the multipath model.
Figure 7 further shows the skyplot coverage rate of four MHMs against the number of continuous days for multipath modeling.The skyplot coverage rate refers to the proportion of grids with at least 30 observations compared to the total number of grids.The MHM_GEC has the highest coverage rate among the four MHMs, surpassing the coverage rate of all the GNSS-specific MHMs (i.e., MHM_G, MHM_E, and MHM_C).However, the coverage rate varies for these MHMs as the number of days increases.The MHM_G coverage rate shows a slow increase with the number of days.This is because the GPS satellite orbit repeats almost every day, so increasing the number of days only marginally improves the coverage.The MHM_C coverage rate increases within 7 days, but then almost stops growing.This is because most of BDS-3's MEOs have an orbit repeat period of 7 days.The coverage rate of MHM_E continues to increase with the number of days until it reaches its orbit repeat period of 10 days.MHM_GEC exhibits a similar growth pattern to MHM_E and MHM_C concerning the number of days.However, MHM_GEC achieves a coverage rate of about 60% with only 3 days of data, while MHM_C requires 7 days.Similarly, MHM_GEC achieves a coverage rate of about 64% with only 4 days of data, while MHM_E takes 10 days to achieve a slightly lower coverage rate of 63%.These results demonstrate that the interoperable MHM (e.g., MHM_GEC) offers higher coverage rate compared to the GNSS-specific MHMs, and can be expected to achieve better compensation effects in a shorter time.

Calibrate carrier-phase residuals on overlap-frequency signals
To assess the efficacy of the four MHMs and ASF model in correcting multipath errors in the observation domain, we constructed these models using 1-10 days of data immediately before the target three days (i.e., days 300, 310, and 320).For brevity without confusion, we hereafter use "n-day" to denote an MHM constructed with n (n = 1, 2 ⋯ 10) days of data.Subsequently, we applied these models separately to GPS-only, Galileo-only, and BDS-3-only PPP solutions to derive carrier-phase residuals with multipath corrections.The reduction rate of residuals can then be calculated by comparing the RMS of residuals before and after the corrections, which serves as a metric of the efficacy of each model in correcting multipath errors in the observation domain.The MHM_G, MHM_GEC, and ASF models fluctuate smoothly with increasing number of days, as indicated by the gray, green, and purple curves in the bottom panels of Fig. 8, because GPS satellites repeat their orbits in space almost every day.Their mean reduction rates of residuals at 12 stations are 24%, 25%, and 21%, respectively.Notably, stations SAS2 and KOS1 experience reduction rates of less than 10%, primarily due to the minimal impact of multipath at these two stations, as explained in the Discussion section.The different designs of satellite orbits result in mean reduction rates of less than 10% and less than 5% for 1-day MHM_E and 1-day MHM_C at the 12 stations, respectively.Nevertheless, as the modeling days increase, their mean reduction rates gradually increase to about 18%.It's worth mentioning that for MHM_C, reduction rates rapidly increase from 1-to 7-day MHM, and reach a plateau from 7-to 10-day MHM_C.This phenomenon can be understood based on the implementation of the 7-day orbit repeat time of BDS-3 MEOs.
Unlike GPS, the Galileo carrier-phase residuals at all 12 stations when corrected with MHM_G exhibit a relatively low reduction rate, as depicted by the light-colored blocks in the left panel of Fig. 9.The mean residual reduction rate at the 12 stations is below 5% corrected with the 1-to 10-day MHM_Gs.Conversely, when corrected with 1-to 10-day MHM_Es, the residual reduction rates display a persistent enhancement at the 12 stations.This is depicted by the gradual transition in color from preblue to dark blue in the second column of panels positioned on the left, proceeding from the left to the right.The corresponding mean reduction rates rapidly increase from 6 to 30%, as indicated by the blue line in the bottom panel.Both MHM_C and MHM_GEC display notable enhancements, with MHM_GEC demonstrating a faster and significantly superior reduction rate.Specifically, the 1-day MHM_GEC exhibits a 12.6% reduction rate which is comparable to the 13.6% reduction rate of the 4-day MHM_E, and surpasses the reduction rate of 5-day MHM_C by 1 percentage point.Moreover, the 6-day MHM_GEC exhibits a reduction rate of 30%, which is merely 0.7 percentage points less than the reduction rate of the 10-day MHM_E and exceeds that of the 10-day ASF by 2 percentage points.The reduction rates of the 7-to 10-day MHM_GECs seem to stabilize at around 31%, exceeding those of other models.Overall, it is evident that the MHM_GEC can achieve faster and more effective enhancements in mitigating multipath effects in comparison to both GNSS-specific MHMs and the ASF.
Like Galileo, MHM_GEC significantly outperforms other models in reducing the BDS-3 multipath effect, especially when a few days of observation are available to construct the MHM, as indicated in Fig. 10.Specifically, the 1-day MHM_GEC exhibits an 11% reduction rate, which surpasses other 1-day MHMs by 5-8 percentage points.The reduction rate increases sharply from the 1-day to 5-day MHM_GEC, with MHM_GEC generally outperforming other MHMs by 8-21 percentage points.It is noteworthy that the 5-day MHM_GEC has a reduction rate of 26%, which is equivalent to that of the 10-day MHM_C, and surpasses that of the 10-day MHM_G, 10-day MHM_E, and 10-day ASF by 22, 3, and 0.5 percentage points, respectively.Additionally, the reduction rates of the 6-to 10-day MHM_GECs increase to around 28.5%, which is 2.5-24.5 percentage points higher than the other 10-day MHMs.Based on these results, it is evident that the multipath mitigation performance of MHM generally improves with the increasing number of modeling days, while the efficiency of the interoperable MHM_GEC is higher.The model can correct multipath effects, which is constructed with only a few days and the correction effect is comparable to that of a GNSS-specific MHM or ASF constructed with one satellite orbit repeat period (10 days for Galileo and 7 days for BDS-3).
The 5-day MHM_GEC demonstrates its excellent performance in mitigating the BDS-3 multipath effect, comparable to or even surpassing the 10-day GNSSspecific MHMs.To further illustrate the specific details of the 5-day MHM_GEC in mitigating the BDS-3 multipath effect, we analyzed a representative time series of BDS-3 C36 satellite carrier-phase residuals before and after correction using various 5-day MHMs, as shown in Fig. 11.The uncorrected carrier-phase residuals exhibit significant multipath variations, such as the peaks and valleys between 20:30 and 22:00 (satellite elevation of 7-15°).After applying various MHMs, the multipath is mitigated to a certain degree.The improvement of MHM_GEC is most significant, reducing the carrierphase residual RMS from 0.44 to 0.30 cm and eliminating most of the multipath fluctuations, as shown with the green dots in the bottom panel.Its carrier-phase residual RMS reduction rate is 31.8%,which is 25, 16, and 9 percentage points higher than that of MHM_G, MHM_E, and MHM_C, respectively.

Positioning improvement after multipath mitigation
To evaluate the performance of positioning using four MHMs and the ASF model, we constructed these models using 1-10 days of data.Subsequently, we applied these models to GPS-only, Galileo-only, and BDS-3only kinematic PPP.The resulting position estimates were then compared with the daily solutions at various stations to calculate the RMS errors.The improvement rate of positioning precision was then calculated by comparing the RMS errors of the PPP solutions with and without model corrections.Since the Galileo constellation is not yet complete, we only considered the periods with more than 5 visible Galileo satellites.The mean improvement rate of positioning precision for the 10 stations (excluding SAS2 and KOS1, which will be Fig.7 Mean skyplot coverage rate (%) of the four MHMs at 12 stations in three time segments against the number of days used.Note that the coverage rate refers to the percentage of grids with a minimum of 30 observations, relative to the total number of grids.The grid resolution is 1° × 1°F ig. 8 Reduction rates of GPS carrier-phase residuals over all three target days corrected by the five models for each station.The color blocks represent the reduction rate of the residual RMS for all GPS satellites in the three target days.The line in the bottom panel represents the mean residual reduction rate over all 12 stations Fig. 9 Reduction rates of Galileo carrier-phase residuals over all three target days corrected by the five models for 12 stations.The color blocks represent the reduction rate of the residual RMS for all Galileo satellites in the three target days.The line in the bottom panel represents the mean residual reduction rate over all 12 stations.The residual reduction rate is shown only for the 10-day ASF due to the 10-day orbital repeat time for Galileo discussed in the next section) on the three target days was calculated to evaluate the performance of each model in the coordinate domain.
Figure 12 illustrates the mean improvement rates of GPS-only (left), Galileo-only (middle), and BDS-3only (right) positioning precision with the calibration of different n-day MHMs and ASF.In the left panels, regarding GPS positioning precision, the 1-day MHM_GEC exhibits significant improvement rates of 38.2%, 42.9%, and 25.5% in the east, north and up components, respectively.However, as the number of days increased, the corresponding improvement rates show only slight variations within 2%, 4%, and 2%.Similarly, the improvement rates of MHM_G and ASF also show slight fluctuations over the days, but their rates are slightly lower than that of MHM_GEC in the east and up components.In contrast, the improvement rates of MHM_E and MHM_C increase significantly as increasing the number of modeling days.The average improvement rate of three-directional positioning precision rises from about 5% for the 1-day model to about 30% for the 10-day model.
In terms of Galileo positioning precision, MHM_GEC outperforms other models in the improvement rate, as depicted in the middle panels of Fig. 12. Specifically, from 1-day to 10-day MHMs, MHM_GEC exhibits the highest improvement rate, followed by MHM_E, MHM_C, and MHM_G.The improvement rate of these MHM models generally increases with increasing number of modeling days, though the growth rate may slow down or reach a plateau.When the number of modeling days reaches 10, the ASF becomes effective with the improvement rates of 39.6%, 42.3%, and 27.8% in the east, north, and up components, respectively, slightly lower than those of the 10-day MHM_GEC and the 10-day MHM_E.Notably, the 6-day MHM_GEC exhibits improvement rates of 35.7%, 40.2%, and 28.1% in the east, north, and up components, respectively, which are very close to those of the 10-day ASF and the 10-day MHM_E.If the 7-day MHM_GEC is utilized, the improvement rates in three components are 38.9%,43.5%, and 29.6%, which are close to the 10-day MHM_E's rates of 39.2%, 43.5%, and 29.9%, and close to the 10-day ASF's rates of 39.6%, 42.3%, and 27.8%.
Similarly, in terms of BDS-3 positioning precision improvement rate, the MHM_GEC outperforms other models, especially when the number of days used to construct MHMs are less than seven.Once the number of days reaches seven, the BDS-3 MEOs repeat period, the ASF becomes effective with the improvement rates of 33.6%, 40.3%, and 19.8% in the east, north, and up components, respectively.These rates are close to the 7-day MHM_C's rates of 35%, 38.2%, and 20.1%, but lower than the 7-day MHM_GEC's rates of 37.5%, 42.5%, and 21.6%.Furthermore, the 5-day MHM_GEC already achieves higher improvement rates of 36.1%, 41.5%, and 21.6% in the east, north, and up components, which are 0.6%, 2.4%, and 1.1% higher than those of 10-day MHM_C and 2.5%, 1.2%, and 1.8% higher than those of 10-day ASF.These results indicate that MHM_GEC is highly efficient in mitigating multipath effects compared to GNSS-specific MHM or ASF models.
Figure 13 displays the reduction rate of horizontal RMS for GPS, Galileo, and BDS-3 kinematic positioning at 10 stations over three target days using the MHMs constructed with five immediately preceding days of data.Across all constellations, the MHM_GEC outperforms other MHMs, reducing the RMS by around 40% in the horizontal component.Notably, the MHM_GEC and the MHM_G exhibit similar reduction rates of about 40% for GPS horizontal positioning error on the three target days, while the MHM_E and the MHM_C show significantly lower mean reduction rates of less than 15%.For Galileo, the mean reduction rates for the MHM_GEC, the MHM_E, the MHM_C, and the MHM_G are 40%, 22%, 18%, and 5%, respectively.For BDS-3, these rates change to 40%, 17%, 14%, and 7%.Despite some variation in reduction rates, the positioning results from all 10 stations demonstrate the advantage of the MHM_GEC over other methods.Overall, the MHM_GEC, which combines overlap-frequency data from GNSS constellations, achieves an impressive improvement of 40% in positioning accuracy even with a short modeling time of 5 days, surpassing the improvement obtained by GNSS-specific MHMs for Galileo and BDS-3.
To further compare the specifics of each MHM in terms of reducing multipath to improve positioning precision, we analyzed a segment that exhibited clear multipath effects from the BDS-3 kinematic PPP solution.Figure 14 shows the time series of BDS-3 kinematic positioning errors in the east, north, and up components before and after the correction of the MHMs constructed with 5 days of data.Without correction, the PPP solution shows significant variations due to multipath errors, but these variations are mitigated to different extents after applying the various MHMs.The most significant improvement is achieved with the MHM_GEC, which reduces the RMS errors in the east, north, and up components from 0.53 cm, 0.8 cm, and 1.0 cm to 0.34 cm, 0.46 cm, and 0.81 cm, respectively.Comparing the top and bottom panels, it is evident that the significant multipath fluctuations in the east and north components are mostly eliminated after the MHM_GEC correction applied.Although the noise in the up component remains relatively high, some large peaks and valleys are visibly reduced, such as the valley between 11:00 and 12:00.On the other hand, the improvement achieved with the MHM_C, the MHM_E, and the MHM_G is not as significant as the MHM_GEC and decreases progressively.Interestingly, the MHM_E reduces the BDS-3 RMS errors with reduction rates of 23%, 30%, and 15% in the east, north, and up components, respectively, which are only 3, 0, and 5 percentage points less than those of the MHM_C.
Figure 15 shows the differences in the mean Power Spectral Density (PSD) of the BDS-3 kinematic positioning errors using different 5-day MHMs at 10 stations.For the periods exceeding 2000s (low-frequency range), the MHM_GEC exhibits the most significant improvement, reducing the positioning error PSD by about 7 dB, 7 dB, and 5 dB in the east, north, and up components, respectively.In the intermediate frequency band (e.g., 2000s to 10 s), PSDs of different MHM solutions slightly increase compared to the solution without MHM correction, such as the MHM_C and the MHM_E.The occurrence of this phenomenon primarily stems from the fact that the multipath correction for each grid is determined by averaging all multipath errors within that grid (Dong et al., 2016).Consequently, it can cause the discrepancies between neighboring grids, which are referred to as "grid effects."This leads to correction jumps when satellites cross over into different grids.Moreover, the duration taken for satellites to pass through a 1° × 1° grid varies from mere seconds to several thousand seconds, thus introducing corresponding jump correction signals in this frequency band.The MHM_GEC excels in correcting multipath, and as more satellite tracks are incorporated into the grid, the variation in correction between adjacent grids Fig. 12 The mean improvement rate (%) of GPS (left), Galileo (middle), and BDS-3 (right) kinematic positioning precision in the east, north, and up components at ten stations on days 300, 310 and 320 after MHM corrections.The horizontal axes represent the number of days used for modeling diminishes, thereby producing only marginal increases in the PSD within this frequency band.Conversely, the MHM_G with its limited satellite track coverage for BDS-3 exhibits a negligible increase and even coincides with the non-corrected solution.In the context of the high-frequency band below 10 s, it is worth noting that all MHMs yield the identical results to the uncorrected solution.This is due to the implementation of a low-pass filter during the experiment, with a cutoff-frequency of 10 s.In conclusion, considering power spectral density, the interoperable MHM exhibits a superior multipath mitigation in comparison to the GNSS-specific MHMs.

Discussion
In Sect.Calibrate carrier-phase residuals on overlap-frequency signals, we observed that stations SAS2 and KOS1 have a smaller multipath correction compared to other stations.Refer to Fig. 16 for better understanding, which shows the hemispherical skyplot of the carrier-phase multipath correction for MHM_GEC at 12 stations.The skyplots at most stations exhibit noticeable color changes, except for stations SAS2 and KOS1.The total RMS of multipath correction for these stations, with elevation ranging from 7° to 90°, is found to be as high as 0.058 cycles.Even at elevations greater than 15° (i.e.15°-90°), there are still blueto-red multipath changes that result in the corresponding multipath correction RMSs of more than 0.02 cycles.In contrast, in the SAS2 and KOS1 multipath skyplots, almost all areas with the elevations greater than 15° are in green, corresponding to carrier-phase multipath corrections within 0.05 cycles, almost at the noise level of the carrier-phase.Only in the low elevation range, i.e. 7°-15°, are there some red and blue shifts.The total RMS multipath correction of these two stations is about 0.02 cycles, which is only 0.3-0.8 of that of the other stations, and the RMS multipath correction for the elevations of 15°-90° is only 0.01 cycles, which is about 0.3-0.5 of that of the other stations.Thus, compared to other stations, the multipath effects at stations SAS2 and KOS1 is weaker.As a result, the improvement in multipath is not significant, regardless of whether the MHM or the ASF is used.

Conclusion
To achieve high-precision GNSS positioning, it is important to mitigate multipath errors.In this study, we compared the positioning precision of GPS, Galileo, and BDS-3 PPP using legacy-frequency signals (GPS L1/L2, Galileo E1/E5a, and BDS-3 B1I/B3I) and overlap-frequency signals (GPS L1/L5, Galileo E1/E5a and BDS-3 B1C/B2a).The results indicate that the positioning precision of GPS and Galileo using legacy-and overlapfrequency signals are almost the same.For BDS-3, the overlap-frequency signals provide a better positioning precision than the legacy-frequency signals, with an improvement of about 18% in the horizontal component and 17% in the vertical component.This demonstrates that using overlap-frequency signals for positioning is feasible and can achieve comparable or better positioning precision than legacy-frequency signals.
Furthermore, we demonstrated the feasibility of creating an interoperable MHM (MHM_GEC) for the overlap-frequency signals of the GPS, Galileo, and BDS-3 satellites.Our evaluation, which compared MHM_GEC's performance with GNSS-specific MHMs and the ASF  and up components at station TIT2 from 8:00 to 12:30 on day 313 of 2021 when using the MHMs constructed with 5 days of data.Note that solutions are shifted vertically by about 6 cm from each other in each direction method, revealed that the MHM_GEC provided the best multipath mitigation performance.The mean reduction of applying the MHM_GEC to GPS, Galileo, and BDS-3 carrier-phase residuals is 25%, 31%, and 28.5%, respectively, which is 3-25 percentage points larger than that of the MHM_G, the MHM_E, and the MHM _C.In particular, we found that the 5-day and 6-day MHM_GECs were comparable to the 10-day MHM_C, the MHM_E and the ASF in reducing BDS-3 and Galileo carrier-phase residuals, and that the 5-to 6-day MHM_GECs resulted in a significant improvement of 40% in positioning precision, comparable to or even exceeding the improvement achieved by the 10-day MHM_E, MHM_C and ASF.
The interoperable MHM_GEC integrates overlapping-frequency observations from different constellations including GPS, Galileo, and BDS-3, which can greatly increase the number of observations to the extent that corrections quickly and uniformly cover the skyplot, resulting in improved spatial resolution, modeling efficiency, and better correction performance in mitigating multipath effects.This means that the interoperable MHM_GEC doesn't need to span a full orbit repeat period as traditional constellation-specific MHM and ASF methods.It can still achieve comparable or even superior multipath mitigation performance.Consequently, this provides an efficient way to model multipath effects for Galileo and BDS-3, enabling rapid response to the changes of station surroundings.

Fig. 1 Fig. 2
Fig. 1 Flowchart for modeling and correcting multipath in MHM and ASF 'G' , 'E' and 'C' stands for GPS, Galileo and BDS-3

Fig. 3
Fig. 3 Distribution of the positioning differences of PPP using legacy-frequency and overlap-frequency signals for all 12 stations over 31 days.The top-left corner of each panel shows the mean and standard deviation (STD) of the positioning differences for all 12 stations over the 31 days Figure 8 presents a comparison of the mean reduction rates of GPS carrier-phase residuals over all three target days corrected by the five models for each station.The color

Fig. 4
Fig. 4 GPS, Galileo, and BDS-3 satellites track in a particular 1° × 1° grid (i.e.48°-49° in elevation and 60°-61° in azimuth) for the overlap-frequency signals at station BAUT from day 310 to day 319 of 2021.The color of the point represents the magnitude of the multipath error

Fig. 6
Fig. 6 The satellite tracks for the overlap-frequency signals at station BAUT on day 319 of 2021.Carrier-phase multipath corrections are color-coded

Fig. 10
Fig. 10 Reduction rates of BDS-3 carrier-phase residuals over all three target days corrected by the five models for 12 stations.The color blocks represent the reduction rate of the residual RMS for all BDS-3 satellites in the three target days.The line in the bottom panel represents the mean residual reduction rate over all 12 stations.The residual reduction rate is shown only for the 7-day ASF due to the 7-day orbit repeat time of BDS-3 MEOs

Fig. 13
Fig. 13 The reduction rate of RMS in horizontal component for GPS, Galileo, and BDS-3 kinematic positioning at 10 stations on days 300, 310, and 320 using the MHMs constructed with five immediately preceding days of data.Each point shows horizontal RMS reduction of a model at a station.Color represents the type of MHMs used.The horizontal line represents the mean RMS reduction of 10 stations

Fig. 14
Fig.14BDS-3 kinematic positioning errors in the east, north, and up components at station TIT2 from 8:00 to 12:30 on day 313 of 2021 when using the MHMs constructed with 5 days of data.Note that solutions are shifted vertically by about 6 cm from each other in each direction

Table 1
Basic information of the stationsSignal and corresponding frequency are for all stations this established model was utilized to generate multipath corrections for dynamic PPP on the target day.During dynamic PPP, station coordinates were estimated as white-noise-like parameters.For the details refer to Table3.The station coordinates are got from the weekly solution coordinate file provided by IGS.All data were calculated using PRIDE PPP-AR, a dedicated open-source software package released by the PRIDE Group

Table 2
Four MHMs established using different constellations in this study

Table 3
Data Processing Strategies and Models

Parameters Strategy or models 1 Hz clock estimation Static PPP Kinematic PPP
W = 1 , e> 30°; W = 4 sine 2 , e < 30° where W is the weight scaling factor and e is the elevation(Geng et al.

Table 4
RMSs of BDS-3 PPP errors using overlap-frequency and legacy-frequency signals

Table 5
RMSs of GPS PPP errors using overlap-frequency and legacy-frequency signals