Correction Model of BDS-3 Satellite Pseudorange Multipath Delays and Its Impact on Single-Frequency Precise Point Positioning

Coal Industry Engineering Research Center of Mining Area Environmental and Disaster Cooperative Monitoring, Anhui University of Science and Technology, Huainan 232001, China School of Environment and Surveying Engineering, Suzhou University, Suzhou 234000, China School of Earth and Environment, Anhui University of Science and Technology, Huainan 232001, China School of Geomatics, Anhui University of Science and Technology, Huainan 232001, China Guangzhou Hi-Target Navigation Tech Co., Ltd., Guangzhou 511400, China


Introduction
BeiDou navigation satellite system (BDS) is a global satellite navigation system independently constructed and operated by China [1]. According to the "three-step" development strategy, the global service has been successfully launched in July 2020 [2][3][4]. As of June 2021, there are more than 50 BDS satellites in orbit, including BDS-2 and BDS-3 satellites, BDS with GPS and Galileo [8]. In addition to providing BDS-3 navigation, positioning, and timing services, B2b can also provide services for international search and rescue and global short message communication [9,10].
In the process of PPP using GNSS, many errors need to be corrected, such as antenna phase center correction and antenna phase wind-up. ese errors can be corrected by the existing algorithms [11][12][13]. However, some errors, such as pseudorange multipath delays, have not yet been corrected by an agreed upon method [14,15]. e downlink navigation signals of satellites are transmitted under nonideal conditions, and the surrounding environment has a great influence on the received signals [2]. When the signal is reflected from the ground and surrounding objects into the receiver interior, the obtained observations contain pseudorange multipath delays and carrier phase multipath delays. e multipath delays of carrier phase observations are generally less than 1/4 wavelength, while the multipath delays of pseudorange observations can reach meter-level. erefore, in some precision measurements, such as mining subsidence and building (structural) health monitoring, the correction of pseudorange multipath delays must be considered [15,16]. e generation mechanism of pseudorange multipath delays is complex, and it has a large relationship with the surrounding reflection environment, so it is difficult to be accurately described [17]. However, the research shows that the pseudorange multipath delay is not a pure random error, but systematic to a certain extent. e systematic part can consider building models to correct the observations, which also makes it possible to deal with the pseudorange multipath delays.
At present, the modeling research of pseudorange multipath delays is mainly in time domain and space domain, and the models obtained by them take time and satellite position as independent variables, respectively. Time domain modeling is considering the periodicity of pseudorange multipath delays, so the application of sidereal filtering in pseudorange multipath delays is mainly studied [16,18,19]. Modeling in space domain is mainly considering the correlation between pseudorange multipath delays and satellite elevation, so the model construction of pseudorange multipath delays based on satellite elevation and azimuth has been studied [2,20]. At present, BDS consists of BDS-2 and BDS-3, including GEO, MEO, and IGSO satellites with different orbital periods, which leads to difficulties in time domain modeling. In contrast, space domain modeling does not require so much periodicity of data and is more flexible. In addition, when the receiver is installed in some dynamic body, the reflection environment of the receiver is consistent, and the space domain method can also be used for modeling.
In 2012, Hauschild et al. found the pseudorange multipath delays in BDS-2 satellite signals for the first time [21]. Subsequently, relevant scholars studied the pseudorange multipath delays of BDS satellite but mostly studied the modeling of BDS-2 satellite, and the pseudorange multipath delays of BDS-3 satellite were less studied [22][23][24]. It is generally believed that the pseudorange multipath delays of BDS-3 satellite are less than those of BDS-2 satellite, which has little influence on PPP. erefore, the impact of pseudorange multipath delays of BDS-3 on PPP is rarely analyzed. By analyzing the pseudorange multipath delays of BDS-3, it is found that the pseudorange multipath delays of BDS-3 can reach decimeter-level or even meter-level under different epochs, which has a certain degree of influence on PPP [25,26].
is paper mainly studied the pseudorange multipath delay of BDS-3 satellite, analyzed its characteristics, and modeled it based on the space domain, and the influence of correcting the pseudorange multipath delays of BDS-3 satellite on single-frequency PPP is analyzed. Firstly, the pseudorange multipath delays of BDS-2 and BDS-3 satellites were compared, and their characteristics were analyzed. en, the EN regularization denoising method was proposed and compared with L2-norm regularization denoising. EN regularization adopts both L2-norm regularization and L1norm regularization to maintain the original characteristics of the sequence but also has a certain sparsity [27][28][29][30]. e pseudorange multipath delays after denoising were modeled by QP + AR model, and the modeling effects were analyzed. Finally, the modeling errors were corrected to the observations, and the positioning accuracy of BDS-3 satellite at each frequency was calculated. e structure of the paper is arranged as follows. Section 2 presents the regularization denoising methods and the modeling principle of QP + AR model. Section 3 discusses the pseudorange multipath delays characteristics of BDS satellite, the regularization denoising effect, the modeling accuracy of QP + AR model, and the PPP accuracy of singlefrequency. Finally, some conclusions and next step works are given in Section 4.

BDS Pseudorange Multipath Delays.
In general, the linear combination of the original pseudorange observations and the carrier phase observations is used to construct multipath observations (MP) to evaluate the multipath delays characteristics of pseudorange observations [2,5].
e MP of the s satellite can be expressed as where the subscripts i and j (i, j � 1, 2, 3; i ≠ j) represent different frequencies, P and φ represent pseudorange observations and carrier phase observations, f, λ denote frequency and wavelength, B is the combination of ambiguity term and hardware delay bias, N is the ambiguity term, d is a constant, and ε denotes other errors. is combination of equations (1) and (2) can eliminate ionospheric delay, tropospheric delay, and geometric errors. When there is no cycle slip, B is generally a constant, which can be obtained by taking the average value of MP.

Mathematical Problems in Engineering
Otherwise, segmented processing is required. When the constant term B is removed, the combined observations can be considered to contain only multipath delays and observation noise. In addition, the multipath delays of the carrier phase observation value is at most 1/4 cycle [16], and the impact is all in centimeter-level, which is negligible compared with the pseudorange multipath delays in decimeter-level or meter-level.

Regularization Denoising
Method. e MP obtained in Section 2.1 contains pseudorange multipath delays and observation noise, so the observation noise should be removed before modeling. Common denoising methods include wavelet threshold denoising [31], empirical mode decomposition [32], and regularization denoising [15,16]. In this work, EN regularization denoising method was proposed and compared with L2-norm regularization denoising. References [15,16] confirmed that L2-norm regularization performed better in multipath delays denoising.

L2-Norm Regularization Denoising Method.
As can be seen from the above content, MP contains pseudorange multipath delays and measurement noise, which can be specifically expressed as where ω represents the weight determined based on the elevation ω e � σ 2 c + σ 2 c /sin 2 (E e ), W � [ω 1 , ω 2 , . . . , ω e , . . . , ω n ], σ c is the pseudorange accuracy, E e is the satellite elevation of e epoch, μ is the Lagrange multiplier, and its solution method can be referred to in [15,16,33].
For the convenience of expression, equation (4) can be expressed in matrix form as follows: with To obtain the optimal valuation of m, equation (5) should satisfy the following equation: Hence, the MP sequence after L2-norm regularization denoising can be obtained as follows:

EN Regularization Denoising Method. According to
Occam's Razor principle, simpler models generally have better generalization performance [29]. erefore, the concept of sparsity in machine learning is introduced here, and the sparse model will be simpler. When the L2-norm regularization is used, the solutions will tend to zero, but not equal to zero, while L1-norm regularization has higher penalty intensity, which can make some solutions equal to zero directly. Considering that L1-norm as a regular term can make the model have better sparsity [33], L1-norm is introduced with L2-norm to form EN regular term together. EN regularization combines the advantages of the two methods, but the computation is also relatively increased. Moreover, both L1-norm regular term and L2-norm regular term can compress the model coefficients.
erefore, EN regularization compresses the model coefficients twice, resulting in a large estimation deviation. However, Zou proposed that a scale transformation can be performed on the EN results to solve this problem [27].
Referring to equation (4), the target function of EN regularization can be constructed as follows: where μ 1 and μ 2 are the hyperparameters of L1-norm regular term and L2-norm regular term, respectively, and other parameters are the same as those in equation (4). For the convenience of expression, equation (9) can be expressed in matrix form as follows: According to equation (10), when μ 1 � 1, μ 2 � 0, the target function contains only L1-norm penalty term, and the EN regularization will be changed to L1-norm regularization. When μ 1 � 0, μ 2 � 1, the target function contains only L2-norm penalty term, and the EN regularization will be changed to L2-norm regularization.
Because the gradient of L1-norm is discontinuous, equation (10) cannot be solved directly [33]; it can be transformed as follows: where τ is a decimal slightly greater than zero. en, equation (10) can be expressed as follows: with ξ (n−1) * (n−1) To obtain the optimal valuation of m, equation (12) should satisfy the following equation: Hence, the MP sequence after EN regularization denoising can be obtained as follows: However, due to the large dimension of the matrix in equation (10), the calculation time is increased. rough analysis, it can be seen that the coefficient matrix is an n-order tridiagonal matrix, which can be expressed as follows: with In equation (15), the matrix elements satisfy the following conditions: erefore, equation (14) can also be solved by the omas algorithm with less computation [15].

MP Model Construction.
e MP of BDS satellite has a certain correlation with the elevation, and the MP sequences after noise reduction can be modeled by QP model. e characteristics of each satellite are not completely consistent, so each satellite is modeled separately. e QP model of the i frequency of satellite s can be expressed as follows: where m is calculated from equation (15), E(e) is the elevation of satellite s in the e epoch, ε is the model noise, and a i0 , a i1 , a i2 are the QP model coefficients, which can be solved by the weighted least square method.
Assuming that the observations of i frequency of S satellites are received in a certain period time, and the polynomial coefficients of all satellites are solved simultaneously, the matrix form of the error equation can be obtained as follows: with After modeling with QP model, the residual presents a stable and random characteristic. AR model is used to further model the residual. e mathematical expression of AR model is as follows: where z s i (e) is the residual value of the satellite s of frequency i after QP modeling in epoch e, p is the order of AR model, ϕ s i,t is the AR coefficient of satellite s on the i frequency, and ε is the white noise.
In a certain period of time, ϕ s i,t can be solved by constructing the error observation equation as follows: e fitting accuracy of AR model is greatly affected by the order of the model. Akaike information criterion (AIC) and Bayesian information criterion (BIC) are commonly used to determine the order of AR model [34,35]. Both AIC and BIC introduced penalty terms related to the number of model parameters to reduce the model parameters and avoid overfitting. However, AIC usually fails when the sample size is large, because the likelihood function value is large, which exceeds the influence of model parameters. In BIC, the punishment related to the number of model parameters is increased, which can effectively control the complexity of the model when the sample size is large. BIC is used to determine the model order, which can be expressed as follows [35,36]: where σ 2 is the mean square error of the model with , p is the number of model parameters, and n is the sample size. BIC value determines the quality of model, and the order with the smallest BIC value is the best order. Due to the large sample size and the amount of computation, the search range can be limited when searching for the optimal order. In this work, we set the order from 1 to 50 to find the optimal order.

Experiments and Analyses
In this section, the data processing method is verified by experiments. e observations of 17 iGMAS stations in 10 days (DOY 1-10, 2021) were collected, and all stations can receive BDS-3 and BDS-2 satellite signals. MEAN, RMS, and Range are selected as statistics to measure accuracy. Since there are few observations of B2b frequency, this work does not analyze its MP, nor study its single-frequency PPP. Figure 1 shows the distribution of iGMAS stations used in the experiment.

BDS MP Analysis.
In this part, the MP of BDS satellites is mainly analyzed. WUH1 station (DOY 1) is selected as an example to visually display the MP characteristics of BDS-2 and BDS-3 satellites at different frequencies.
e MP characteristics of each satellite at each frequency are analyzed. In this section, one satellite of each type is selected as an example, which is shown in Figure 2 e MP of BDS-2 satellite changes with the elevation more obviously than that of BDS-3 satellite. e MP of different frequencies of the same satellite has obvious difference. For BDS-2 satellite, the influence of B3I is significantly less than that of B1I. Range and RMS show that the influence of B3I is more stable, and the fluctuation amplitude is smaller than that of B1I. For BDS-3 satellite, the fluctuation amplitude of MP at four frequencies shows the rule of B1C > B1I > B3I > B2a. e MP of B2a is relatively minimum, but it can still reach 0.5 m∼1.0 m at low elevation. B1C has the largest MP fluctuation amplitude. e MP of GEO satellite of BDS-2 and BDS-3 and MEO satellite of BDS-2 shows obvious periodicity at B1I and B3I frequencies. However, MEO and IGSO of BDS-3 did not show obvious periodicity. From the perspective of RMS and Range, the fluctuation amplitude of GEO satellite of BDS-2satellite is lower than that of MEO and IGSO satellite, while that of BDS-3 is opposite.

MP Denoising and Modeling.
e regularization algorithm in Section 2.2 is used to denoise the MP sequence. Take C28 and C40 in Figure 2 as examples; both satellites contain four frequencies. For the convenience of distinction, the original MP sequence is denoted as MP, the L2-norm regularization denoising sequence is denoted as MP-regL2, and the EN regularization denoising sequence is denoted as MP-regEN. Figure 3 shows the comparison of the two satellites before and after the L2-regularization and EN regularization. Figure 4 shows the RMS and Range variations of the two satellites after MP regularization at each frequency. Table 1 shows the RMS statistical results of different satellites at different frequencies after regularization denoising.
It is obvious that the two regularization methods used in the experiment have the ability to denoise. e regularized sequence is still consistent with the volatility of the original MP sequence, but the sequence is relatively more stable. EN regularization makes the sequence sparse to a certain extent, and some small data in the original sequence become zero after EN regularization. In addition, the size of the hyperparameter is closely related to the stability of the sequence. e larger the hyperparameter is, the more stable the sequence is, but the data are prone to distortion. erefore, it is extremely important to determine the size of the hyperparameter reasonably. Table 1 shows that the stability of MP sequences at different frequencies of different satellites is improved after regularization, and the stability after EN regularization is better than that after L2-norm regularization.
In this section, the modeling method in Section 2.3 is used to model MP. Most of the existing results suggest that the MP of BDS-2 MEO and IGSO satellites is significantly correlated with the elevation, while the MP of BDS-3 satellite is weakly correlated with the elevation, which makes it impossible to establish an effective model. In this work, QP + AR model is used to model the MP of BDS-3 satellite and analyze the modeling effect. Due to the small number of GEO satellites and the small range of elevation variation, only MEO and IGSO satellites of BDS-3 are studied in this part. e MP sequences after L2-norm regularization and EN regularization are modeled, respectively, and the model coefficients of each satellite at each frequency are calculated,  Mathematical Problems in Engineering respectively. Figures 5 and 6 show the modeling results for C28 and C40 in Figure 2.
Overall, the results of QP + AR model are basically consistent with the original sequence.
e results of QP model mainly reflect the change of trend terms. In most cases, the impact on BDS-3 is within ±0.2 m, and the contribution to the modeling results is relatively small. By analyzing the coefficients of QP model, it is found that QP fitting parameters are relatively small. Table 2 lists the QP model parameters of C28 and C40. It is found that the coefficients of the QP model are small; especially, the coefficients of a 2 are close to zero. AR model is more suitable for the modeling of stationary random sequences, which can accurately represent the fluctuation of sequences. However, it should be noted that the order of AR model is fixed during the modeling process. e order of AR adaptively selected by different epochs should be studied in the following work.

Single-Frequency PPP Experiment.
In order to further test the modeling effect of Section 3.2, two iGMAS stations (WUH1 and KUN1; DOY 1-10, 2021) are selected for BDS-3 single-frequency PPP experiment. e precision ephemeris and clock bias products are mainly downloaded from MGEX. Table 3 shows some data processing strategies. ree single-frequency PPP schemes are designed. Scheme 1: correct the errors without MP, and then carry out single-frequency PPP according to Table 3. Scheme 2: similar to scheme 1, but also correct MP. Using L2-norm regularization to denoise MP, the QP + AR is established, and the modeling results are corrected to the observation value, and then the single-frequency PPP is carried out according to Table 3. Scheme 3: similar to scheme 2, but EN regularization is used instead of L2-norm regularization.
To intuitively show the positioning accuracy of BDS-3 single-frequency PPP of the three schemes, the positioning    Figures 7 and 8. Tables 4 and 5 show the positioning accuracy statistics of the two stations after convergence. It can be seen that when BDS-3 satellite is used only for single-frequency PPP, the four frequencies in the E and N directions can reach centimeter-level after convergence, but the U direction is basically in the decimeter-level. When the observation environment of the station is consistent, the B2a frequency of BDS-3 performs best in single-frequency PPP.
Compared with Scheme 1, Scheme 2 and Scheme 3 show that the overall positioning accuracy is improved. Increasing MP correction can significantly improve the positioning accuracy of B3I and B2a frequencies, in which the improvement of the E and U directions can reach 10.6%-34.4% and 5.9%-65.7%, and the improvement of the N direction accuracy is mostly less than 10.0%. e positioning accuracy of B1I and B1C frequencies can be improved by 0%∼30.7% and 0.4%∼28.6% in the E and U directions, respectively, while the improvement in the N direction is small, and even the accuracy will be reduced.
Compared with Scheme 2, Scheme 3 performs better in the single-frequency PPP of B3I and B2a, while Scheme 2 performs better in the single-frequency PPP of B1I and B1C. e denoising of EN regularization is greater than that of L2norm regularization, and it can make the sequence show a certain degree of sparsity. erefore, the fluctuation amplitude of MP after denoising is lower than that of L2-norm regularization, which also leads to different QP + AR modeling values of the same epoch. As can be seen from Section 3.1, the MP of B1I and B1C is greater than that of B3I and B2a, so the MP correction of B1I and B1C frequency in Scheme 2 is greater than that in Scheme 3. e MP of B3I and B2a is relatively small, and the fluctuation is more stable, and in this case, the correction effect of Scheme 3 is better.
Convergence time is an important indicator of PPP performance. Generally, the positioning deviation in E, N, and U directions is less than 0.1 m and maintained for a certain time as the convergence condition [39]. In this work, BDS-3 single-frequency observation data is used for PPP,  but the positioning accuracy is slightly lower, especially in the U direction. Based on the positioning results of many tests and considering the large positioning deviation in the U direction, the convergence condition in this work is defined as: the positioning deviation in the E and N directions is better than 0.1 m, in the U direction is better than 0.6 m, and  the positioning deviation in the next 20 epochs is within the limit. Figure 9 shows the convergence time of WUH1 and KUN1 stations in DOY 1.
As can be seen from Figure 9, the convergence time of each frequency is quite different, among which B1I frequency has the longest convergence time, B3I frequency has the second, and B2a frequency has the shortest convergence time. In terms of convergence time and accuracy, the accuracy of the new signals transmitted by BDS-3 is better than that of the old signals, and the quality of B2a frequency observation data is relatively better. For the three schemes, the convergence time after MP correction is shortened to a certain extent, but it is not obvious. However, compared with Scheme 1 and Scheme 3, the convergence time of B1I frequency of KUN1 station is significantly shorter in Scheme 2. It is found that the positioning deviation of Scheme 1 and Scheme 3 has also changed smoothly after 7h, but the U direction is slightly greater than 0.6, so it is not counted as  e convergence time of B1I frequency of WUH1 station is longer, which is also because the positioning deviation of the U direction is not less than 0.6 m. However, the positioning deviation of the three schemes also tends to be stable in about 6h. is is also the case for the positioning deviation of other frequencies. Due to the low accuracy of single-frequency PPP, how to set the convergence threshold rationally needs further study.

Conclusions
Pseudorange multipath delay is one of the main errors that affect the PPP, and an accurate correction model can improve the positioning accuracy. is work studied the influence of BDS-3 satellite MP on single-frequency PPP. Firstly, the characteristics of each frequency MP of BDS satellite are analyzed, and then the MP sequence is denoised by L2-norm regularization and EN regularization, respectively. e QP + AR model is constructed for the denoised sequence, and the modeling value is corrected to BDS-3 single-frequency PPP. e following conclusions can be drawn from the experiments: (1) e MP of BDS-3 satellite can reach the decimeterlevel. Although it is generally smaller than that of BDS-2 satellite, the influence cannot be ignored. Among the four frequencies, the MP of B1C is the largest, and that of B2a is the smallest. However, when the elevation is relatively low, the MP of BDS-3 can still reach 0.5 m∼1.0 m, which also indicates that the MP of BDS-3 is correlated with the elevation to a certain extent.
(2) EN regularization is a combination of L1-norm regularization and L2-norm regularization. Compared with the L2-norm regularization denoising, MP sequences have certain sparsity and higher stability. However, it should be noted that the size of  the regularization parameter has a great influence on the denoising effect, and its size needs to be determined reasonably.
(3) e QP + AR is constructed by using the denoised MP and corrected to the pseudorange observations, which can improve the accuracy of single-frequency PPP. e accuracy improvement of B2a and B3I is more obvious than that of B1I and B1C, and the improvement in the E and U directions is significantly greater than that in the N direction. e accuracy of the E and U directions can be improved by 10.6%∼34.4% and 5.9%∼65.7%, and the N direction can be improved by less than 10.0%, but in some cases, the accuracy can even be reduced. (4) When B3I and B2a use EN regularization denoising sequence to construct QP + AR model to correct MP, their single-frequency PPP accuracy is better than that of L2-norm regularization denoising modeling to correct MP, while B1I and B1C use L2-norm regularization denoising modeling to correct MP with higher single-frequency PPP accuracy. (5) When BDS-3 single-frequency PPP converges, the position deviation in E and N directions can reach centimeter-level, and the U direction remains in demerit-level. e signal quality of B2a and B1C is better than that of B1I and B3I in terms of convergence time and accuracy.
In general, it is feasible to denoise MP and correct it by modeling, and the accuracy of BDS-3 single-frequency PPP can be improved to a certain extent by correcting MP. However, the number of stations used in this work is not large, which cannot be universal as the model parameters in reference [20]. In addition, the order of AR model is fixed when calculating model coefficients. In the subsequent research, it can be considered to change the order adaptively in different epochs. Although more computation is required, better modeling accuracy will be achieved. When the order of AR model is too large, we can also consider adding regularization to prevent overfitting. In addition, the positioning accuracy of single-frequency PPP is slightly lower, especially in the U direction, so how to reasonably determine the convergence threshold also needs to be further studied.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare no conflicts of interest.