The Extraction and Characterization of Pseudorange Multipath Based on BDS-3 Multi-Frequency Observations

Global Navigation Satellite System (GNSS) observations are subject to various errors during their propagation process. A reasonable correction of these errors can improve the positioning, navigation, and timing (PNT) service capability. The impact of multipaths on pseudorange observations can reach a decimeters or even meters level. However, their mechanism is complex and there is currently no universally accepted high-precision correction model. The correlation between the pseudorange multipaths (MP) of BDS-2 satellites and satellite elevation has been confirmed, while there have been fewer analyses of the MP characteristics for different frequencies of BDS-3 satellites. The broadcasting of multi-frequency observations in BDS-3 should theoretically make the extracted MP more accurate compared to traditional methods. Based on this, in this contribution, a multi-frequency MP extraction algorithm based on the least squares principle is proposed, which can simultaneously eliminate the influence of higher-order ionospheric delay. The analytical expression for only eliminating first-order ionospheric delay is successfully derived. Subsequently, the characteristics of the MPs extracted from different frequency combinations and the impact of combination noise on the extraction accuracy are discussed. The influence of second-order ionospheric delay on the MPs for each frequency under different combination noises, as well as the periodic behavior exhibited in long-term observations of the BDS-3 medium earth orbit (MEO) and inclined geosynchronous orbit (IGSO) satellites, are also analyzed. Finally, the correlations between the MPs of each frequency of BDS satellite and elevation are quantitatively analyzed based on observations from 35 stations. Overall, this work has positive implications for the study of the MP characteristics of BDS-3 and subsequent modeling efforts.


Introduction
GNSS has been widely used in PNT services due to its advantages of high accuracy, high frequency, and all-weather capabilities [1,2]. As the main member of navigation systems, the Beidou Satellite Navigation System (BDS) was independently designed and developed by China and has undergone a "three-step" development strategy. On 31 July 2020, it was officially opened for global service. Compared with other navigation systems, the BDS satellite constellation is composed of three types of orbiting satellites: geostationary orbit (GEO), IGSO, and MEO. Both the IGSO and MEO satellites of BDS-3 can transmit observations of five frequency signals [3,4]. Despite the increased complexity in process-

Methods
Carrier phase and pseudorange observations are the basic observations broadcasted by each navigation system and their linearized observation equations can be expressed as [35]: Φ i = ρ s r + cdt r − cdt s + T r − k i I 1,1 − ν i I 1,2 + λ i N i + ε Φ i P i = ρ s r + cdt r − cdt s + T r + k i I 1,1 + 2ν i I 1,2 + MP i + ε P i (1) where the subscript i (i = 1, 2, . . . , n) represents different frequencies; Φ and P represent the carrier phase and pseudorange observations; ρ s r denotes the geometric distance between the satellite and the receiver; dt r and dt s represent the receiver and satellite clock errors, c is the speed of light; T r represents the tropospheric delay; k i = f 2 1 / f 2 i and ν i = f 3 1 / f 3 i represent the first-and second-order ionospheric delay amplification factors of frequency i, respectively, I 1,1 and I 1,2 denote the first-and second-order ionospheric delay of the f 1 frequency, λ i and N i represent the wavelength and ambiguity of frequency i, MP i is the pseudorange multipath, and ε Φ i and ε P i denote other errors of the carrier phase and pseudorange observations.
Usually, MP i can be extracted from dual-frequency observations and its expression is [26,31,36]: where i and j denote different frequencies, B i,j is the combination of the ambiguity and hardware delay, and ξ represents other errors. When the carrier phase observation does not include cycle slip, B i,j can be treated as a constant, that is, after removing B i,j from the sequence extracted from Equation (2), an MP sequence with certain noise can be obtained [37]. The magnitude of the noise can have a certain impact on the modeling and correction of an MP. As shown in Equation (2), the extraction of MPs at different frequencies only uses the carrier phase observations of two frequencies. However, with the increase in the frequencies broadcasted by various navigation systems, introducing additional frequency carrier phase observations should be able to obtain MP sequences with less noise. For BDS-3 MEO and IGSO satellites, five signal frequencies can now be broadcasted. If fivefrequency carrier phase observations are used to extract an MP at a certain frequency, it can be represented as: where η 0 and η j represent the combination coefficients of the pseudorange and carrier phase observations at frequency j, respectively, while η i denotes the combination coefficients of the carrier phase observations at other frequencies.
From Equation (1), it follows that, for the MP 1 extraction using Equation (3), when only eliminating the influence of I 1,1 , η 0 ∼ η 5 should satisfy: When the effects of I 1,1 and I 1,2 are considered simultaneously, Equation (4) can be combined with Equation (5).
According to Equations (4) and (5), η 1 ∼ η 5 can be obtained by the Lagrange multiplier method. In fact, the solution of Equation (4) should introduce additional parameters, which increases the computational workload. In this work, starting from the observation Equations of the pseudorange and carrier phase, a method for calculating the multi-frequency MP combination coefficients based on the least squares principle is proposed, directly according to the characteristics of the combined parameters.
Equation (4) shows that the coefficients of the MP combination need to satisfy the condition that the combined geometric distance term and ionospheric delay term are both equal to zero. Therefore, when considering only the geometric distance term, I 1,1 , and MP, the observation equations for the pseudorange and carrier phase can be simplified as follows:  where ρ denotes the sum of the geometric distances of the all frequency-independent terms from the satellite to the receiver. Assuming that the covariance matrix of the observation is Q, the least squares solution of the estimated parameters X can be obtained from Equation (6), which is expressed as: Equation (7) can also be seen as a combination of the pseudorange observations and carrier phase observations at different frequencies to extract ρ, I 1,1 , and MP 1 , where the combination coefficients are: In order to satisfy the condition of Equation (4), the transformation matrix At this point, the combined ρ and I 1,1 can be kept at 0, while the coefficient of MP 1 remains at 1, that is: According to Equation (9), the combination coefficient of MP 1 can be expressed as: The optimal coefficient of the MP combination for a certain frequency can be obtained through Equation (10). Assuming that the carrier phase observations of n frequencies are used, the MP combination coefficient can be expressed as: Equations (11) and (12) provide analytical expressions for the optimal coefficients of multi-frequency MP combinations. When extracting MPs using different frequency observations, it is only necessary to directly bring in the k i corresponding to different frequencies. However, it is important to note that k i is not fixed and constant. Although the same carrier phase observations are used, when the frequency of extracting the MP is different, k i also needs to be updated again.
If the influence of I 1,2 needs to be eliminated when extracting an MP, Equation (6) can be rephrased as: By updating the transformation matrix R(R 4 = 0 0 0 1 ) accordingly, the MP combination coefficients eliminating I 1,2 can be quickly obtained. Taking triple-frequency observations as an example, the combination coefficient when simultaneously eliminating I 1,1 and I 1,2 can be expressed as: Taking the B1I frequency as an example, Table 1 shows the combination coefficients and carrier phase observation noise amplification factors (Ω = η 1 2 + η 2 2 + · · · + η n 2 ) for extracting MPs using different frequency combinations. The DF, TF, QF, and FF represent dual-frequency, triple-frequency, quad-frequency, and five-frequency, respectively. Figure 1 shows the combination noise of the B1C/B1I/B3I/B2a frequency when extracting MPs through different frequency combinations.  From Table 1 and Figure 1, it can be seen that the Ω varies significantly under different frequency combinations. For example, the Ω of B1I/B1C can be 25-35 times that of other dual-frequency combinations and about 45 times that of quad-and five-frequency combinations. Related studies have pointed out that, when the Ω is large, it will amplify the observation error implied by the carrier phase observations [20,34]. Therefore, combinations of Ω similar to B1I/B1C should be careful not to be used for extracting MPs. The commonly used B1I/B3I combination has a combined noise of 6.245. From Figure 1, it can be seen that the Ω is greater than most combinations. In theory, when the Ω is smaller, the extracted MP sequence contains less noise. Therefore, combinations such as triple-frequency or quad-frequency should be able to extract MP sequences with smaller error fluctuations.
Additionally, when both the I1,1 and I1,2 are eliminated simultaneously, the Ω significantly increases and the optimal combination is different from when only I1,1 is considered. For example, the combination of B1I/B1C/B2a for B1I is relatively optimal for the triple-frequency combination. However, when both I1,1 and I1,2 are simultaneously considered, this combination actually has the largest combination noise. Through the analysis of the coefficients of the four frequency combinations, it is found that the contribution of the B3I frequency to the combination is relatively small when considering only I1,1, while it has the largest contribution when eliminating both I1,1 and I1,2 simultaneously. Another point to note is that, when extracting MPs with the same frequency of different frequencies, the combination coefficients are different. Although the expression of the coefficients is the same, due to the change in the first frequency, the corresponding k2 and k3 will also From Table 1 and Figure 1, it can be seen that the Ω varies significantly under different frequency combinations. For example, the Ω of B1I/B1C can be 25-35 times that of other dual-frequency combinations and about 45 times that of quad-and five-frequency combinations. Related studies have pointed out that, when the Ω is large, it will amplify the observation error implied by the carrier phase observations [20,34]. Therefore, combinations of Ω similar to B1I/B1C should be careful not to be used for extracting MPs. The commonly used B1I/B3I combination has a combined noise of 6.245. From Figure 1, it can be seen that the Ω is greater than most combinations. In theory, when the Ω is smaller, the extracted MP sequence contains less noise. Therefore, combinations such as triple-frequency or quad-frequency should be able to extract MP sequences with smaller error fluctuations.
Additionally, when both the I 1,1 and I 1,2 are eliminated simultaneously, the Ω significantly increases and the optimal combination is different from when only I 1,1 is considered. For example, the combination of B1I/B1C/B2a for B1I is relatively optimal for the triplefrequency combination. However, when both I 1,1 and I 1,2 are simultaneously considered, this combination actually has the largest combination noise. Through the analysis of the coefficients of the four frequency combinations, it is found that the contribution of the B3I frequency to the combination is relatively small when considering only I 1,1 , while it has the largest contribution when eliminating both I 1,1 and I 1,2 simultaneously. Another point to note is that, when extracting MPs with the same frequency of different frequencies, the combination coefficients are different. Although the expression of the coefficients is the same, due to the change in the first frequency, the corresponding k 2 and k 3 will also change. For instance, when using the combination of B1C/B1I/B2a to extract the MP of each frequency, the corresponding combination noises are Ω B1C = 3.593, Ω B1I = 3.621, and Ω B2a = 4.820, respectively.

Experiments and Analyses
In this section, the correctness of the theories proposed in this work and the characteristics of the MPs for each frequency of BDS-3 will be discussed through experiments. B3I, and B2a frequencies of BDS-3. The sampling interval for the observations at each station was 30 s. Figure 2 shows the distribution of these stations. Additionally, a set of dynamic vehicle observations (CAR1 station) was also selected to further validate the theories proposed in this work. The data acquisition scheme can be found in reference [38]. Standard deviation (STD), mean, and range were selected as statistics to measure the accuracy.

Experiments and Analyses
In this section, the correctness of the theories proposed in this work and the characteristics of the MPs for each frequency of BDS-3 will be discussed through experiments. The observations of 35 stations (18 MGEX stations and 17 iGMAS stations) over 7 days (DOY 1-7, 2021) were selected, and all stations could receive observations of the B1I, B1C, B3I, and B2a frequencies of BDS-3. The sampling interval for the observations at each station was 30 s. Figure 2 shows the distribution of these stations. Additionally, a set of dynamic vehicle observations (CAR1 station) was also selected to further validate the theories proposed in this work. The data acquisition scheme can be found in reference [38]. Standard deviation (STD), mean, and range were selected as statistics to measure the accuracy.

Comparison of Characteristics of MP Extracted with Different Frequency Combinations
In order to compare the differences in the extracted MP sequences when different frequency combinations were used, taking the B1I and B1C frequencies as examples, four modes with varying combination noise from dual-frequency to quad-frequency were selected to extract the MPs for each frequency. The specific combination frequencies and coefficients are shown in Table 2. Although there were differences in the combined noise of the four modes, each mode listed was the optimal combination coefficient calculated for the corresponding frequency. This section mainly verifies the difference between the different frequency combinations on the obtained MP sequences, without considering the difference between the optimal combination and non-optimal combination among each frequency. In fact, the experimental results can still reflect the difference between them, as the combination noise of the non-optimal combinations also increased correspondingly, and the comparison results should be similar to the results.

Comparison of Characteristics of MP Extracted with Different Frequency Combinations
In order to compare the differences in the extracted MP sequences when different frequency combinations were used, taking the B1I and B1C frequencies as examples, four modes with varying combination noise from dual-frequency to quad-frequency were selected to extract the MPs for each frequency. The specific combination frequencies and coefficients are shown in Table 2. Although there were differences in the combined noise of the four modes, each mode listed was the optimal combination coefficient calculated for the corresponding frequency. This section mainly verifies the difference between the different frequency combinations on the obtained MP sequences, without considering the difference between the optimal combination and non-optimal combination among each frequency. In fact, the experimental results can still reflect the difference between them, as the combination noise of the non-optimal combinations also increased correspondingly, and the comparison results should be similar to the results.
Taking the observation of the KUN1 station as an example, Figure 3 presents the MP sequences extracted from the C28 and C40 satellites using different frequencies. In each sub-figure, the upper portion represents the extracted MP sequences for the B1I and B1C frequencies, while the lower portion shows the first-order difference between each mode and the QF mode. Table 3 provides the accuracy statistics for the MP sequences extracted using different combination modes, as well as the improvement relative to DF1.
By comparing the four modes, it can be observed that the MP sequences extracted using the triple-frequency and quad-frequency combinations exhibited a similar trend to those extracted using the traditional dual-frequency method, but with smaller fluctuations. This indicates that the coefficient calculated using Equation (11) was correct and optimal. The results also show that if the Ω was small, the fluctuation of the extracted MP sequence was smaller. By comparing this with the QF mode, it can be observed that, when the difference of the Ω was relatively small, the first-order difference between the two sequences generally fluctuated within 3 cm. However, the difference between the QF mode and the DF1 mode still fluctuated at the meter level. This indicates that, when the Ω was small, the extracted sequences were essentially consistent with each other. Additionally, Table 3 also shows that the improvement in the accuracy, relative to the DF1 mode, was consistent across the different combination modes. Therefore, when certain frequency observations were missing, other combination modes could be adaptively selected to extract the MP. Figure 4 presents histograms of the frequency distributions for the MP sequences extracted using different modes. It can be seen that the sequences from the DF2, TF, and QF modes are relatively more concentrated around zero. However, due to the favorable observation environment at the KUN1 station, the representation in the histograms of the frequency distribution is not significant.  Taking the observation of the KUN1 station as an example, Figure 3 presents the MP sequences extracted from the C28 and C40 satellites using different frequencies. In each sub-figure, the upper portion represents the extracted MP sequences for the B1I and B1C frequencies, while the lower portion shows the first-order difference between each mode and the QF mode. Table 3 provides the accuracy statistics for the MP sequences extracted using different combination modes, as well as the improvement relative to DF1.    mode, was consistent across the different combination modes. Therefore, when ce frequency observations were missing, other combination modes could be adaptivel lected to extract the MP. Figure 4 presents histograms of the frequency distribution the MP sequences extracted using different modes. It can be seen that the sequences the DF2, TF, and QF modes are relatively more concentrated around zero. However to the favorable observation environment at the KUN1 station, the representation i histograms of the frequency distribution is not significant. To further validate the theories proposed in this article and compare the charac tics of the MP sequences extracted from the different modes, the same frequency co nations and coefficients shown in Table 2 were used to process the dynamic observat Similar to Figure 3, Figure 5 presents the results and difference sequences for sate C23 and C38, while Figure 6 displays histograms of the frequency distribution for th sequences extracted using the different modes. Table 4 provides the accuracy statistic the two satellites.  To further validate the theories proposed in this article and compare the characteristics of the MP sequences extracted from the different modes, the same frequency combinations and coefficients shown in Table 2 were used to process the dynamic observations. Similar to Figure 3, Figure 5 presents the results and difference sequences for satellites C23 and C38, while Figure 6 displays histograms of the frequency distribution for the MP sequences extracted using the different modes. Table 4 provides the accuracy statistics for the two satellites.     It can be seen that, due to the significant influence of the surrounding environment during the collection of the dynamic observations, the selected multi-frequency combinations in this work could effectively reduce the fluctuation level and magnitude of the MP sequences. Figure 5 provides a more intuitive representation of the numerical variation range of the MPs extracted using different combinations compared to the DF1 mode. Although the improvement in the accuracy achieved using the DF2, TF, and QF modes for the dynamic observations was not significantly different from that of the DF1 mode, this accuracy improvement was more pronounced compared to the KUN1 station. The B1I and B1C frequencies of the two satellites at the KUN1 station showed accuracy improvements ranging from 38-42% to 13-15%, respectively, while the CAR1 station exhibited even higher improvements, reaching over 30-39% for the B1I frequency and 43-50% for the B1C frequency. The comparison between the two stations highlights the increased importance of using optimal multi-frequency combinations in dynamic experiments.

Analysis of the Influence of Second-Order Ionospheric Delay
In Section 3.1, the influence of the four different frequency combination modes on the MP extraction was compared using the B1I and B1C frequencies as examples. However, the coefficients used were the optimal solutions only when I 1,1 was eliminated, and the influence of the simultaneous elimination of I 1,2 on the extracted MPs was not considered. From the derivation process of the optimal coefficients in Section 2, it can be found that, when only I 1,1 was considered, the dual-frequency combination could only have a unique solution, while the triple-frequency, quad-frequency, and five-frequency combinations could obtain the optimal solution with the minimum combination noise due to the presence of redundant observations. Similarly, when considering both I 1,1 and I 1,2 , Equations (4) and (5) need at least three frequency observations to achieve a unique solution, and the quad-frequency and five-frequency combinations could obtain the optimal solution. Based on the findings in Section 3.1, when the combination noise was similar, the extracted MP sequences were also very similar. In this section, a further analysis is conducted to examine the impact of ionospheric delay on these MP sequences. Taking the triple-frequency and quad-frequency combinations as examples, the influence of second-order ionospheric delay on the MP sequences extracted using different frequency combinations is analyzed. Table 5 provides the selected combination modes and their corresponding coefficients for the B1C, B3I, and B2a frequencies. The combination modes for the B1I frequency can be found in Table 1, where B1I/B1C/B3I, B1I/B1C/B2a, B1I/B3I/B2a, and B1I/B1C/B3I/B2a are represented as B1ITF1, B1ITF2, B1ITF3, and B1IQF1, respectively. Theoretically, regardless of whether the I 1,2 was eliminated or not, the coefficients of the pseudorange observations involved in the combination were all 1, which meant that the truth values contained in the extracted sequences were consistent. However, due to the different combination coefficients of the two methods, the errors contained in the extracted sequence also differed. In order to analyze the impact of eliminating I 1,2 on the extraction of the MP sequences, the optimal coefficients of the same mode in Table 5 were used to extract the sequences for eliminating I 1,1 and I 1,1 /I 1,2 , and a difference was made between the two sets of sequences. The difference between the two sets of sequences could be considered as the influence of I 1,2 under that mode. Taking the continuous 7-day observations from the WUH2 station as an example, Figures 7 and 8, respectively, show the difference sequences for C28 (MEO satellite) and C38 (IGSO satellite) at various frequencies. B1CTF1 in the figures represents the first-order difference sequence for the B1CTF1 mode, and the meanings of the other sequences are similar. Table 6 provides statistical information on the accuracy of all the BDS-3 MEO and IGSO satellites at the experimental stations for the 7 days under each mode. be considered as the influence of I1,2 under that mode. Taking the continuous 7-day ob vations from the WUH2 station as an example, Figures 7 and 8, respectively, show difference sequences for C28 (MEO satellite) and C38 (IGSO satellite) at various frequ cies. B1CTF1 in the figures represents the first-order difference sequence for the B1C mode, and the meanings of the other sequences are similar. Table 6 provides statis information on the accuracy of all the BDS-3 MEO and IGSO satellites at the experime stations for the 7 days under each mode.  be considered as the influence of I1,2 under that mode. Taking the continuous 7-day ob vations from the WUH2 station as an example, Figures 7 and 8, respectively, show difference sequences for C28 (MEO satellite) and C38 (IGSO satellite) at various freq cies. B1CTF1 in the figures represents the first-order difference sequence for the B1C mode, and the meanings of the other sequences are similar. Table 6 provides statis information on the accuracy of all the BDS-3 MEO and IGSO satellites at the experime stations for the 7 days under each mode.    Figures 7 and 8, as well as Table 6, it can be observed that:  Based on Figures 7 and 8, as well as Table 6, it can be observed that: (1) There was a clear correlation between the magnitude of the fluctuations in the difference sequences and the ∆Ω. The ∆Ω for the B1ITF1, B1ITF2, B1CTF1, B1CTF2, B3ITF1, and B2aTF1 modes were all above 140. The STD and range values for each sequence were also greater than 0.24 m and 1.5 m, respectively, which were 3-5 times larger than those for the other modes. Among them, the B2aTF1 mode had the largest ∆Ω, which was 245.626. Its STD and range were 0.401 m and 2.478 m, respectively, which were also significantly higher than those of the other modes. With the exception of the MP extraction mode for the B2a frequency, the fluctuation range of the differenced sequences for other modes was less than 0.5 m and the STD was basically around 0.05-0.07 m.
(2) According to reference [32,33], when only I 1,1 is considered, the anti-multipath performance of the B2a frequency is better than that of other BDS-3 frequencies. However, the experiments showed that, when eliminating I 1,2 simultaneously, the MPs of the B2a frequency extracted using different modes exhibited larger fluctuations. Table 6 also indicates that, except for modes such as B2aTF1 with significant noise, the STDs for the other three modes were approximately 2-3 times larger than those for the other frequency modes. From the statistical results, it can be seen that the averages of the STDs and ranges of the B2aTF2, B2aTF3, and B2aQF1 modes were 0.113 m and 0.726 m, respectively, which were about 1.8, 1.7, and 2.2 times and 1.6, 1.7, and 1.9 times greater than those of the B1I, B1C, and B3I frequency combinations. The anti-multipath performance of the B3I frequency was also relatively superior, and the combination with a smaller ∆Ω corresponded to a lower STD and range compared to those of the other frequencies. This characteristic was contrary to that of the B2a frequency, and further research is needed in the future.
(3) By comparing the differenced MP sequences of the C28 and C38 satellites at various frequencies, it could be found that, regardless of the magnitude of the fluctuations in the differenced sequences, the differenced sequences of the C38 satellite at each frequency exhibited a clear periodicity. However, the periodicity of the C28 satellite was not prominent. Compared with other IGSO and MEO satellites, they all showed the same periodicity pattern. Additionally, when comparing the four frequencies, it could be observed that their periodic patterns were consistent.
In this section, the analysis of the MP difference sequences under the different modes was also conducted using the dynamic observations from Section 3.1. Figures 9 and 10 provide the MP difference sequences extracted with the consideration of I 1,2 for various frequencies of the C23 and C38 satellites at the CAR1 station. Table 7 presents the accuracy statistical information of all the MEO and IGSO satellites at the CAR1 station using different modes.
According to the experimental results, it can be concluded that the fluctuation magnitude of the differenced MP sequences of the BDS-3 IGSO satellites at the CAR1 station was smaller than that of the MEO satellites. The STD and range for each frequency were approximately half of those for the MEO satellites. In addition, similar to the WUH2 station, the fluctuation amplitude of each sequence was positively correlated with the ∆Ω. The fluctuation magnitude of the differenced sequences for the B2a frequency in various modes remained higher than that of the B1I, B1C, and B3I frequencies and the STD and range were also around 1.5-2.0 times larger. The accuracy of the B3I frequency was also better than that of the other frequencies. Due to the observation time being only 2 h, all the sequences fluctuated around the zero value, and the mean values were in the millimeter or sub-millimeter level. The IGSO satellites did not exhibit periodicity. Compared with the WUH2 station, although the observation at the CAR1 station was more affected by the environment, overall, the magnitude of the second-order ionospheric delay remained consistent. For modes with a high ∆Ω, the STD and range of the CAR1 station were slightly higher than those of the WUH2 station. However, when the ∆Ω was small, the STD and range were actually better than that of the WUH2 station. their periodic patterns were consistent.
In this section, the analysis of the MP difference sequences under the different mode was also conducted using the dynamic observations from Section 3.1. Figures 9 and 1 provide the MP difference sequences extracted with the consideration of I1,2 for variou frequencies of the C23 and C38 satellites at the CAR1 station. Table 7 presents the accurac statistical information of all the MEO and IGSO satellites at the CAR1 station using diffe ent modes. pattern. Additionally, when comparing the four frequencies, it could be observed th their periodic patterns were consistent. In this section, the analysis of the MP difference sequences under the different mod was also conducted using the dynamic observations from Section 3.1. Figures 9 and provide the MP difference sequences extracted with the consideration of I1,2 for vario frequencies of the C23 and C38 satellites at the CAR1 station. Table 7 presents the accura statistical information of all the MEO and IGSO satellites at the CAR1 station using diffe ent modes.

Analysis of the Correlation between MP and Elevation of BDS-3 Satellite
Based on the previous experiments, it is known that the impact of the BDS-3 satellite's MP on the pseudorange observations can reach the order of decimeters or even meters, and it needs to be corrected during precise positioning. Currently, there have been numerous studies discussing MP correction methods for BDS-2 satellites, but relatively fewer studies have been conducted for BDS-3 satellites [30,39]. Referring to the relevant research on BDS-2, this section focuses on analyzing the relationship between the MPs and satellite elevations for BDS-3 MEO and IGSO satellites. It should be noted that this analysis only examined their correlation and did not investigate the modeling algorithms, as that would require other processing methods, which will be addressed in future research. For the analysis of this correlation, the MP sequences extracted using the QF mode in Section 3.1 were utilized. Furthermore, in order to avoid MPs being applied to specific satellites, uniform processing was applied to a specific type of satellite at each station to reflect its statistical characteristics. As of February 2023, the distribution of BDS-3 MEO and IGSO satellites can be found in Table 8. The satellites marked with bold numbers indicate their current status as on-orbit testing and do not participate in statistics.  Taking IGSO satellites as an example, assuming that a station observes a total of n IGSO satellites during a certain period and that the effective dataset of the f i frequency for each satellite is X n f i = {X 1 , X 2 , · · · , X n−1 , X n }, the dataset composed of the satellite elevation (ele) and MP for all satellites and all epochs is Q IGSO BDS−2 = ele(X n By using this method, a dataset consisting of elevation and MP can be obtained for each type of satellite and frequency, with both the elevation and MP containing 900 elements. Considering that different stations receive satellites at different elevations and that the maximum and minimum elevation for receiving a particular satellite may vary, there may be cases in the dataset where the MP is empty for certain elevation angles. In such cases, the maximum elevation in Q(d k ) can be adaptively set based on the specific satellite reception conditions at each station. Figures 11 and 12 provide examples of the statistical results of the MPs and satellite elevations for MAYG and OWMG stations, respectively.
The analysis of the statistical results for all the stations in Figure 2 reveals that the MPs of the MEO and IGSO satellites of BDS-3 exhibited similar systematic biases related to elevation, similar to BDS-2 satellites. Both of them showed more significant changes in MEO satellite systematic biases with elevation. The correlation strength between the MP and elevation varied significantly among different types of satellites and frequencies.
The anti-multipath performance of BDS-3 at various frequencies was better than that of BDS-2. The MP of BDS-3 fluctuated around −0.5 m-0.5 m, while BDS-2 mainly fluctuated around −1.0 m-1.0 m. Comparing the different frequencies, it can be seen that the MP fluctuation amplitude of the B2a and B3I frequencies was smaller, while that of the B1C and B1I frequencies was larger. Additionally, the magnitude of the MP fluctuations varied across different satellite orbits for the same frequency. Overall, for BDS-3 satellites, the MP magnitude of the MEO and IGSO satellites was roughly equivalent, while for BDS-2 satellites, the magnitude of the MEO satellites was greater than that of the IGSO.
responding to the fi frequency when the height angle is dk' = dk − 0.05°. By using this method, a dataset consisting of elevation and MP can be obtained for each type of satellite and frequency, with both the elevation and MP containing 900 elements. Considering that different stations receive satellites at different elevations and that the maximum and minimum elevation for receiving a particular satellite may vary, there may be cases in the dataset where the MP is empty for certain elevation angles. In such cases, the maximum elevation in ( ) k Q d can be adaptively set based on the specific satellite reception conditions at each station. Figures 11 and 12 is formed. The average of ( ) k MP d is taken as the MP corresponding to the fi frequency when the height angle is dk' = dk − 0.05°. By using this method, a dataset consisting of elevation and MP can be obtained for each type of satellite and frequency, with both the elevation and MP containing 900 elements. Considering that different stations receive satellites at different elevations and that the maximum and minimum elevation for receiving a particular satellite may vary, there may be cases in the dataset where the MP is empty for certain elevation angles. In such cases, the maximum elevation in ( ) k Q d can be adaptively set based on the specific satellite reception conditions at each station. Figures 11 and 12   In order to quantify the correlation between the MP and elevation for each frequency, the correlation coefficient was calculated for different frequencies at each station. The calculation equation is [40]: where X and Y are the two vectors for the correlation to be determined, X and Y are the average values corresponding to each vector, and n is the length of the vector. The value of r xy is between −1 and 1, with r xy > 0 indicating a positive correlation and r xy < 0 indicating a negative correlation. Generally, the absolute of r indicates the strength of the correlation. Typically, |r|≥ 0.8 indicates a very strong correlation, 0.6 ≤|r|< 0.8 indicates a strong correlation, 0.4 ≤|r|< 0.6 indicates a moderate correlation, 0.2 ≤|r|< 0.4 indicates a weak correlation, and |r|< 0.2 is uncorrelated [41]. Table 9 shows the average correlation coefficient of each experimental station. From Table 9, it can be observed that, for GEO satellites, the absolute values of the correlation coefficients for each frequency ranged from 0.06 to 0.13, indicating a clear lack of correlation. For IGSO and MEO satellites, the absolute values of the correlation coefficients for each frequency ranged from 0.42 to 0.56 and 0.63 to 0.80, respectively, showing moderate and strong correlations. Additionally, the correlation between the MEO and IGSO satellites of BDS-2 was stronger at the B1I and B3I frequencies than that of BDS-3. Among the four frequencies of BDS-3, the B2a frequency exhibited the strongest correlation.

Summary and Conclusions
The broadcasting of the multi-frequency observations in various navigation systems presents new opportunities for high-precision PNT services. However, the pseudorange and carrier phase observations of each frequency are subject to various interferences during signal propagation. As one of the errors affecting the pseudorange observations, MPs have also been studied by many scholars. On the basis of existing research, the MP characteristics of BDS-3 satellites at various frequencies were studied based on multi-frequency observations in this work. Through the analysis of static and dynamic observations, the main conclusions can be drawn as follows: (1) Starting from the pseudorange and carrier phase observation equations, a direct multi-frequency MP combination coefficient calculation method based on the least squares principle was proposed. This method was simple in its calculation and could effectively utilize the observation information of each frequency. On this basis, a detailed formula for calculating the MP combination coefficients for eliminating I 1,1 and I 1,2 was derived, and a unified analytical expression was summarized when only eliminating I 1,1 . Compared with the traditional dual-frequency MP combination extraction results, the trend of the MP sequence extracted by the TF and QF combination was basically consistent with it, and the fluctuation amplitude was smaller, which confirmed the correctness of the proposed method.
(2) The accuracy of the MP sequences extracted from the four modes with different Ω was compared by using static and dynamic observations. The results showed that the magnitude of Ω directly affected the fluctuation amplitude of the extracted MP sequences. When there was a significant difference in Ω, the difference in the extracted MP sequences could reach the decimeter level. When the difference of Ω was small, the difference in the extracted MP sequences was generally within 3 cm. Compared to the DF1 mode, the modes with a smaller Ω showed an average improvement of over 25% in the accuracy of the extracted MP sequences, with an even higher accuracy improvement for dynamic observation. Additionally, the selected combinations were the optimal combination coefficients for each frequency combination, which also indirectly reflected the comparison of the accuracy between the optimal and non-optimal combinations in the MP extraction.
(3) Taking the triple-frequency and quad-frequency combinations as examples, the influence of I 1,2 on the extraction of the MPs for the different frequencies and modes was analyzed by taking the first-order difference with the MP sequence obtained by only eliminating I 1,1 . The fluctuation amplitude of the differenced MP sequence also showed a clear correlation with ∆Ω, with the STD and range differing by 3-5 times. The B2a and B3I frequencies exhibited better anti-multipath performances. However, under the same level of combination noise, the STD and range of the differenced MP sequence of the B2a frequency were both more than 1.5 times larger than those of the other three frequencies, while the STD and range of the B3I frequency were the smallest.
(4) The impact of I 1,2 on the IGSO and MEO satellites of BDS-3 was different. In different modes, the first-order difference sequences of the four frequencies of the IGSO satellites always exhibited a significant periodicity, while the periodicity of the MEO satellites was not significant. In addition, for the static observation, the magnitude of the differenced MP sequences between the two types of satellites was basically consistent, but for the dynamically collected observation, the STD and range of the IGSO satellites were about half those of the MEO satellites.
(5) Based on the observations from all the experimental stations, the correlation between the MPs of the IGSO and MEO satellites and the satellite elevation was analyzed. The MP of the MEO satellites showed a strong correlation with the elevation, while the correlation for the IGSO satellites was of moderate strength. The correlation between the MEO and IGSO satellites of BDS-2 at the B1I and B3I frequencies was stronger than that of BDS-3. Among the four frequencies of BDS-3, the B2a frequency exhibited the strongest correlation. Overall, the anti-multipath performance of BDS-3 frequencies was superior to that of BDS-2.
This work focused on the method of extracting MPs using multi-frequency observations and analyzed the characteristics of MPs under different scenarios. However, in practical positioning processes, it is necessary to consider how to correct such errors to improve the positioning accuracy. This involves addressing issues such as the extraction of the true information from MPs and the construction of high-precision correction models. These aspects are important considerations for multi-frequency and multi-system precise positioning and will be one of the future research directions.