# Improved Mitigation Method for the Multipath Delays of BDS-3 Code Observations with the Aid of a Sparse Modeling Algorithm

Authors: Chao Hu; Qianxin Wang; Alberto Hernandez Moraleda
Journal: Journal of Sensors
 (2021)

---

## Abstract

Global navigation satellite systems are essential for positioning, navigation, and timing services. The quality and reliability of satellite observations determine the system performance, especially in the case of the newly launched global BDS-3 service. However, analyses of multipath delays in BDS-3 satellite observations suggest that there are appreciable errors at different frequencies. Improvement of the accuracy and precision of positioning, navigation, and timing services provided by BDS-3 requires the mitigation of multipath delays of the satellite observations. This paper models the multipath delays of BDS-3 observations using a least-squares combined autoregressive method. Furthermore, a sparse modeling algorithm is proposed to obtain a multipath delay series using total variation and elastic net terms for denoising and eliminating the effect of limited original observations. The estimated coefficients of multipath delays are then set as prior information to correct the next-arc code observations, where the square-root information filter is used in the coefficient estimation. Moreover, four groups of experiments are conducted to analyze the results of modeling the BDS-3 multipath delay using the proposed methods, with single-frequency precise point positioning (PPP) and real-time PPP solutions being selected to test the correction of multipath delays in BDS-3 code observations. The residuals of iGMAS and MGEX station coordinates indicate improvements in eastward, northward, and upward directions of at least 4.1%, 9.6%, and 1.2%, respectively, for the frequency B1I; 6.6%, 5.3%, and 0.2%, respectively, for B3I, 12.5%, 14.3%, and 3.8%, respectively, for B1C; and 5.9%, 7.4%, and 18.1%, respectively, for B2a relative to the use of the traditional method in BDS-3 single-frequency PPP. Furthermore, the real-time double-frequency PPP is optimized by at least 10% forB1I+B3I and B1C+B2a. An improved result was obtained with the proposed strategy in a standard point positioning experiment. The proposed multipath delay mitigation method is therefore effective in improving BDS-3 satellite code observations.


---

## Body

## 1. Introduction

Global navigation satellite systems (GNSSs) provide all-weather and continuous services to global users in the area of high-precision spatial information. China has developed the BeiDou Navigation Satellite System (BDS) in a three-step strategy involving a verification system (BDS-1), regional system (BDS-2), and global system (BDS-3) [1, 2]. BDS-3 was officially announced to be operational in 2020, representing a shift in the BDS from regional to global services [3]. At present, more than 40 BDS satellites in the orbit provide positioning, navigation, and timing (PNT) services (http://www.csno-tarc.cn). BDS-3 will provide the majority of spatial and temporal references for services in the future. Therefore, the high performance of BDS-3 services, especially in terms of the precision and accuracy of observations, is essential for the development of BDS [4, 5].The GNSS observation quality relates to the ambiguity resolution [6], bias parameters [7, 8], and correction models [9]. The parameter estimation of the satellite orbit and clock offsets is affected by errors in the satellite-ground observations, such as multipath delays. Meanwhile, the precision of precise point positioning (PPP) is directly limited by the quality of satellite observations [10, 11]. In terms of quality control, studies on BDS observations have focused on the C/N0 ratio [12, 13], corrections for noise and multipath delay [14], and precision analyses of the phase and code [15]. Moreover, for BDS-3 multifrequency observations, the characteristics of the observation signals [16], qualities of new observations [17], and biases of data fusion [3] have typically been discussed. In summary, the quality of BDS-3 satellite observations is better than that of traditional BDS-2 satellite observations [16]. However, the results of experiments on BDS-3 code multipath delays suggest that there is a trend term in the code observations [2]. Therefore, to mitigate BDS-3 code multipath delays, methods of modeling BDS-3 multipath delays need to be studied for high-precision applications.In the processing of GNSS observation multipath delays, there are four aspects to high-precision estimation and modeling methods: (1) time domain filtering algorithms, such as the wavelet algorithm and Vondrak algorithm, (2) modeling using the repeating periods of satellites, such as modeling using the sidereal filtering algorithm [18], (3) compensation algorithms in the spatial domain, which provide a multipath hemispherical map [14] and allow spherical harmonic analysis [19], and (4) machine learning methods, such as Tikhonov regularization and sparse modeling [20, 21]. The literature indicates that all the methods mentioned above are useful for the mitigation of BDS code multipath delays. Notably, new signals and their modulation and multiplexing techniques are used for BDS-3 satellites. However, widely used strategies for BDS code multipath delays focus on BDS-2 observations, for example, a two-step modeling method, which first models the multipath delay and then the code bias [2], cannot be applied to BDS-3. Owing to differences in the signal characteristics, signal numbers, signal quality, and satellite spatial distributions between BDS-2 and BDS-3, BDS-3 observations should be carefully investigated when using multiple frequencies. The present study divides the code multipath delays of BDS-3 observations (at four frequencies) into systematic and random parts for modeling and mitigation, respectively.It is difficult to directly construct the multipath delay with a theoretical model owing to the effects of the complex environment and receiver front end. Additionally, the estimated multipath delays are mainly contaminated by noise, which cannot be modeled and mitigated for using an empirical algorithm [20]. Additionally, the Tikhonov regularization algorithm has been used to denoise a series of BDS clock offsets [22] and model GNSS multipath errors [20], with the results of experiments suggesting that the regularization algorithm obtains a better PPP solution. Moreover, Kalman filtering [23], spectral analysis [24], and the particle filter [25] have been proposed to denoise a series of GNSS multipath delays. Meanwhile, to improve the accuracy of models, machine learning methods, such as sparse modeling [21] and the kernel trick [26], have been presented recently for processing GNSS observations, where the least absolute shrinkage selection operator (Lasso), elastic net, and fast iterative shrinkage threshold algorithm (FISTA) are used to output the estimated model coefficients. However, two important issues, namely, the length of the multipath delay series and the condition of noise in the observations, are ignored in widely used sparse modeling strategies [27]. In general, elastic net is an improved version of the Lasso algorithm to conduct the sparse modelling, which is more stable for the extreme correlations among various predictors as the combination of L1 norm and L2 norm regularization terms. Thus, elastic net is set in the process of minimizing the loss function, which ensures the sparsity of the model and increases the stability of the model solution. Additionally, total variation uses the nonsmooth norm to replace the traditional L2 norm regularization, in which the local structure information can be captured to further eliminate the influence of noise on the reconstruction of a complex network structure and to increase the robustness of the model. In consideration of the characteristic of BDS-3 code observations and the limited sample points, sparse modeling combined with the elastic net term EN and the total variation term TV [27] is investigated for improving the robustness of BDS-3 multipath modeling in the present paper.The aim of this paper is at optimizing the model of BDS-3 code multipath delay and at mitigating the effects of the delay on rapid PPP solutions. A new strategy for modeling the BDS-3 code multipath delay is proposed in Section2, where sparse modeling combined with EN and TV terms is first presented, an integrated estimation of model coefficients is then derived in detail, and the correction model of next-arc observations is finally designed using prior information. Section 3 tests single-frequency PPP, real-time PPP, and standard point positioning (SPP) using the proposed strategy. Conclusions and perspectives are presented in Section 4.


## 2. Materials and Methods

The service performance of a satellite system directly reflects observation multipath delays as a main error in GNSS PPP. To reduce the effects of noise and limited observations, this section discusses the high-precision modeling of BDS-3 code multipath delays, for which the systematic and random parts are separately considered. Meanwhile, a strategy of correcting for the code multipath delay of the next-arc observation is also studied.


### 2.1. Modeling BDS-3 Code Multipath Delays

To overcome the effects of the observation noise and series length in modeling BDS-3 multipath delays, sparse modeling from the field of machine learning [27] is adopted to optimize and process the series of BDS-3 multipath errors (raw code multipath delays for each satellite with continuous ambiguity calculated by equations (1) and (2)). In this research, based on the continuous tracking of satellites, the multipath delays of different BDS-3 code observations are first estimated using a combination of two frequencies as [2].
(1)MPistk=Pistk−mijs⋅λi⋅ϕistk+mijs−1⋅λj⋅ϕjstk,(2)mijs=fi2+fj2fi2−fj2,where f and λ are the frequency and wavelength of the signals, respectively, the subscripts i and j indicate different frequencies, s represents the satellite (BDS-3 in this paper), φ and P are the observations of the phase and code, respectively, and tk is the epoch time.Using the estimated code multipath delays, a series function is written for theith BDS-3 satellite as
(3)yit1yit2⋮⋮⋮yitM=φi1t1φi2t1⋯φiNt1φi1t2φi2t2⋯φiNt2⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮φi1tMφi2tM⋯φiNtMxi1xi2⋮⋮⋮xiN+ςi1ςi2⋮⋮⋮ςiM,where y denotes the multipath delays of different epochs or elevations (estimated values), ɸ denotes the independent variable of multipath models, which can be constructed using the epoch, elevation, and azimuth, x denotes the model coefficients, ζ denotes the residuals of the models, and M and N are the total number of epochs and the length of the series, respectively.Equation (3) is rewritten in matrix form as
(4)Yi=φi⋅Xi+ςi.Then, to improve the robustness of equation (4) against observation noise, EN of the sparse structure and the TV term are introduced into the objective function [28]:
(5)FMi=Yi−φiXi22+β2φiXi22+λφiXi1+γTVφiXi=min,where β/2φiXi22+λφiXi1 denotes the regularization term of EN, TVφiXi is the regularization term of TV, and β, λ, and γ are coefficients of the objective function, which can be set as prior values or selected using the Thomas algorithm [20]. The smoothness strategy of TV regularization is adopted to obtain solutions to equation (5) [29] as follows:
(6)TVφiXi=∑A⋅φiXi,A=−11−11⋯⋯−11−11.We adopt the concept of FISTA to iterate the results of model coefficients using the inequality approximation algorithm and the corresponding interval threshold of convergence (ε) [30]. The final solutions of equation (5) are thus obtained using the iteration condition
(7)FMik+1−FMik≤ε,where k is the number of iterations.


### 2.2. Integrated Modeling of BDS-3 Code Multipath Delays

On the basis of the sparse modeling of BDS-3 code observation multipath delays, the PPP service can be improved theoretically through the precise correction of multipath errors. However, the short-term prediction of multipath delays should be obtained for next-arc or real-time PPP uses, where a more precise model coefficient is needed to construct the BDS-3 code multipath delay. Adopting the integrated estimation strategy, the systematic and random parts of the multipath series are separated and modeled in one solution. Taking thekth BDS-3 satellite as an example, the code multipath delays are expressed as a function of the satellite elevation:
(8)Ykei=Bkei+Rkei+εkei,where ei is the ith elevation of the satellite, Yk denotes the code multipath delays after processing with equations (1)–(7), and Bei and Rei are functions of systematic and random parts, respectively. Traditionally, the systematic part is modeled with a polynomial function, where a second-order function is used to process the BDS code multipath delays [2]. The random part can be expressed as an autoregressive (AR) model. ε represents the model residuals.As mentioned above, the AR model is used to describe the random terms of multipath errors, where the relationship among variables of a random series is expressed as(9)Rkei=κk1Rkei−1+⋯+κkdRkei−1+μkei.Here,κ is the coefficient of the AR model estimated by the least-squares (LS) algorithm, d is the order of the model, and μ is the white noise. In equation (9), the order of the AR model is selected using the Akaike information criterion [22, 29]. To estimate the model coefficients properly and rapidly, we take all coefficients into one solution as
(10)Ykei+1=Kkei+1⋅bk+Rkei+1+υkei,where υ is the model error, Kei=1eiei2, and bk=b0kb1kb2kT. Equation (10) can be rewritten as
(11)Rkei+1=Ykei+1−Kkei+1⋅bk+υkei=κk1⋅Ykei−Kkei⋅bk+⋯+κkd⋅Ykei+1−d−Kkei+1−d⋅bk+υkei.Furthermore, inserting equation (10) into equation (11) yields
(12)Ykei+1−κk1⋅Ykei+⋯+κkd⋅Ykei+1−d+υk′ei=Kkei+1⋅bk−∑h=1dκkh⋅Kkei+1−h⋅bk+υk′ei.Stacking all equations for satellite elevations fromef to et, the model coefficients for the kth satellite are obtained from the integrated estimation as
(13)cd+1Tcd+2T⋮ctT−F⋅Nd+1Nd+2⋮Nt⋅Yk=Aked+1Aked+2⋮Aket−F⋅Md+1Md+2⋮Mt⋅bk,,where ⊗ denotes the Kronecker product, F=It−d⊗κkT, κk=κk1κk2⋯κkdT, It−d is an identity matrix with order t–d, cj is the jth column of the t×t identity matrix, Nj=cj−1Tcj−2T⋯cj−dTT, and Mj=Akej−1Akej−2⋯Akej−dT.The model coefficients are obtained from the estimation of equation (13). Using equation (11), all satellites (e.g., satellite l) are processed in one solution:
(14)C⋅Y′=T⋅b,(15)Y′=Y1Y2⋯Yk⋯YlTb=b1b2⋯bk⋯blC=C1C2⋯Ck⋯ClT=T1T2⋯Tk⋯TlCk=cd+1Tcd+2T⋯ctTT−F⋅Nd+1Nd+2⋯NtTTk=Aked+1Aked+2⋯AketT−F⋅Md+1Md+2⋯MtT.


### 2.3. Multipath Delay Correction of Next-Arc Observations

The estimated coefficients are traditionally introduced into next-arc observations to mitigate the multipath delays [2]. However, considering the orbital period and the surrounding environments of the receivers, it is difficult to obtain a perfect PPP solution under the effect of code multipath delays. For correction of the multipath delays of the next-arc BDS-3 code observations using the estimated coefficients, a new strategy is proposed to further improve the accuracy of multipath mitigation. The estimation of the next-arc multipath delays using the LS algorithm is expressed as
(16)r=Bx+ηP,where x denotes the true values of the model coefficients for the next-arc multipath delay and P is the weight matrix. In addition, to improve the accuracy of the model coefficients, a constraint function based on the prior information x0, x~0=x+ε~0, is used. Here, a weight matrix P~0 is set. ε and η are the model residuals. During the estimation of equation (16), the weight matrix is decomposed as P~0=R~0−1R~0−T and P=R−1R−T.Equation (16) and its constraint function are then rewritten as
(17)z~0=R~0⋅x+v~0,z=Ax+v,where z~0=R~0⋅x~0, z=R⋅v, A=R⋅B, v~0=R~0⋅ε~0, and v=R⋅η.Meanwhile, the Householder conversion for equation (17) is further expressed as
(18)S⋅R~0Ax=S⋅z~0z+S⋅v~0v⇒R^00x=z^0e+v^0ve,where e and ve are the model errors and the corresponding noise, respectively.An objective function is defined as(19)Jx=R~0Ax−z~0z2=R∧0x−z∧02+e2≥e2.In obtaining the minimumJx, the parameter solution reads as
(20)R^0x=z^0.Finally, the solution to equation (20) is
(21)x=R^0−1⋅z^0,where the covariance matrix of the estimated x is P^0=R^0−1⋅R^0−T.According to equations (16)–(21), the code multipath delay of BDS-3 can be further mitigated, especially in rapid applications, with the aid of prior information. Figure 1 is a flowchart of the integrated method of processing BDS-3 code observations proposed in this research.Figure 1
Flowchart of mitigating BDS-3 code observation multipath delays.


## 2.1. Modeling BDS-3 Code Multipath Delays

To overcome the effects of the observation noise and series length in modeling BDS-3 multipath delays, sparse modeling from the field of machine learning [27] is adopted to optimize and process the series of BDS-3 multipath errors (raw code multipath delays for each satellite with continuous ambiguity calculated by equations (1) and (2)). In this research, based on the continuous tracking of satellites, the multipath delays of different BDS-3 code observations are first estimated using a combination of two frequencies as [2].
(1)MPistk=Pistk−mijs⋅λi⋅ϕistk+mijs−1⋅λj⋅ϕjstk,(2)mijs=fi2+fj2fi2−fj2,where f and λ are the frequency and wavelength of the signals, respectively, the subscripts i and j indicate different frequencies, s represents the satellite (BDS-3 in this paper), φ and P are the observations of the phase and code, respectively, and tk is the epoch time.Using the estimated code multipath delays, a series function is written for theith BDS-3 satellite as
(3)yit1yit2⋮⋮⋮yitM=φi1t1φi2t1⋯φiNt1φi1t2φi2t2⋯φiNt2⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮φi1tMφi2tM⋯φiNtMxi1xi2⋮⋮⋮xiN+ςi1ςi2⋮⋮⋮ςiM,where y denotes the multipath delays of different epochs or elevations (estimated values), ɸ denotes the independent variable of multipath models, which can be constructed using the epoch, elevation, and azimuth, x denotes the model coefficients, ζ denotes the residuals of the models, and M and N are the total number of epochs and the length of the series, respectively.Equation (3) is rewritten in matrix form as
(4)Yi=φi⋅Xi+ςi.Then, to improve the robustness of equation (4) against observation noise, EN of the sparse structure and the TV term are introduced into the objective function [28]:
(5)FMi=Yi−φiXi22+β2φiXi22+λφiXi1+γTVφiXi=min,where β/2φiXi22+λφiXi1 denotes the regularization term of EN, TVφiXi is the regularization term of TV, and β, λ, and γ are coefficients of the objective function, which can be set as prior values or selected using the Thomas algorithm [20]. The smoothness strategy of TV regularization is adopted to obtain solutions to equation (5) [29] as follows:
(6)TVφiXi=∑A⋅φiXi,A=−11−11⋯⋯−11−11.We adopt the concept of FISTA to iterate the results of model coefficients using the inequality approximation algorithm and the corresponding interval threshold of convergence (ε) [30]. The final solutions of equation (5) are thus obtained using the iteration condition
(7)FMik+1−FMik≤ε,where k is the number of iterations.


## 2.2. Integrated Modeling of BDS-3 Code Multipath Delays

On the basis of the sparse modeling of BDS-3 code observation multipath delays, the PPP service can be improved theoretically through the precise correction of multipath errors. However, the short-term prediction of multipath delays should be obtained for next-arc or real-time PPP uses, where a more precise model coefficient is needed to construct the BDS-3 code multipath delay. Adopting the integrated estimation strategy, the systematic and random parts of the multipath series are separated and modeled in one solution. Taking thekth BDS-3 satellite as an example, the code multipath delays are expressed as a function of the satellite elevation:
(8)Ykei=Bkei+Rkei+εkei,where ei is the ith elevation of the satellite, Yk denotes the code multipath delays after processing with equations (1)–(7), and Bei and Rei are functions of systematic and random parts, respectively. Traditionally, the systematic part is modeled with a polynomial function, where a second-order function is used to process the BDS code multipath delays [2]. The random part can be expressed as an autoregressive (AR) model. ε represents the model residuals.As mentioned above, the AR model is used to describe the random terms of multipath errors, where the relationship among variables of a random series is expressed as(9)Rkei=κk1Rkei−1+⋯+κkdRkei−1+μkei.Here,κ is the coefficient of the AR model estimated by the least-squares (LS) algorithm, d is the order of the model, and μ is the white noise. In equation (9), the order of the AR model is selected using the Akaike information criterion [22, 29]. To estimate the model coefficients properly and rapidly, we take all coefficients into one solution as
(10)Ykei+1=Kkei+1⋅bk+Rkei+1+υkei,where υ is the model error, Kei=1eiei2, and bk=b0kb1kb2kT. Equation (10) can be rewritten as
(11)Rkei+1=Ykei+1−Kkei+1⋅bk+υkei=κk1⋅Ykei−Kkei⋅bk+⋯+κkd⋅Ykei+1−d−Kkei+1−d⋅bk+υkei.Furthermore, inserting equation (10) into equation (11) yields
(12)Ykei+1−κk1⋅Ykei+⋯+κkd⋅Ykei+1−d+υk′ei=Kkei+1⋅bk−∑h=1dκkh⋅Kkei+1−h⋅bk+υk′ei.Stacking all equations for satellite elevations fromef to et, the model coefficients for the kth satellite are obtained from the integrated estimation as
(13)cd+1Tcd+2T⋮ctT−F⋅Nd+1Nd+2⋮Nt⋅Yk=Aked+1Aked+2⋮Aket−F⋅Md+1Md+2⋮Mt⋅bk,,where ⊗ denotes the Kronecker product, F=It−d⊗κkT, κk=κk1κk2⋯κkdT, It−d is an identity matrix with order t–d, cj is the jth column of the t×t identity matrix, Nj=cj−1Tcj−2T⋯cj−dTT, and Mj=Akej−1Akej−2⋯Akej−dT.The model coefficients are obtained from the estimation of equation (13). Using equation (11), all satellites (e.g., satellite l) are processed in one solution:
(14)C⋅Y′=T⋅b,(15)Y′=Y1Y2⋯Yk⋯YlTb=b1b2⋯bk⋯blC=C1C2⋯Ck⋯ClT=T1T2⋯Tk⋯TlCk=cd+1Tcd+2T⋯ctTT−F⋅Nd+1Nd+2⋯NtTTk=Aked+1Aked+2⋯AketT−F⋅Md+1Md+2⋯MtT.


## 2.3. Multipath Delay Correction of Next-Arc Observations

The estimated coefficients are traditionally introduced into next-arc observations to mitigate the multipath delays [2]. However, considering the orbital period and the surrounding environments of the receivers, it is difficult to obtain a perfect PPP solution under the effect of code multipath delays. For correction of the multipath delays of the next-arc BDS-3 code observations using the estimated coefficients, a new strategy is proposed to further improve the accuracy of multipath mitigation. The estimation of the next-arc multipath delays using the LS algorithm is expressed as
(16)r=Bx+ηP,where x denotes the true values of the model coefficients for the next-arc multipath delay and P is the weight matrix. In addition, to improve the accuracy of the model coefficients, a constraint function based on the prior information x0, x~0=x+ε~0, is used. Here, a weight matrix P~0 is set. ε and η are the model residuals. During the estimation of equation (16), the weight matrix is decomposed as P~0=R~0−1R~0−T and P=R−1R−T.Equation (16) and its constraint function are then rewritten as
(17)z~0=R~0⋅x+v~0,z=Ax+v,where z~0=R~0⋅x~0, z=R⋅v, A=R⋅B, v~0=R~0⋅ε~0, and v=R⋅η.Meanwhile, the Householder conversion for equation (17) is further expressed as
(18)S⋅R~0Ax=S⋅z~0z+S⋅v~0v⇒R^00x=z^0e+v^0ve,where e and ve are the model errors and the corresponding noise, respectively.An objective function is defined as(19)Jx=R~0Ax−z~0z2=R∧0x−z∧02+e2≥e2.In obtaining the minimumJx, the parameter solution reads as
(20)R^0x=z^0.Finally, the solution to equation (20) is
(21)x=R^0−1⋅z^0,where the covariance matrix of the estimated x is P^0=R^0−1⋅R^0−T.According to equations (16)–(21), the code multipath delay of BDS-3 can be further mitigated, especially in rapid applications, with the aid of prior information. Figure 1 is a flowchart of the integrated method of processing BDS-3 code observations proposed in this research.Figure 1
Flowchart of mitigating BDS-3 code observation multipath delays.


## 3. Experiments and Results

To verify the proposed strategy of correcting for BDS-3 observation multipath delays, BDS-3 observations of the Multi-GNSS Experiment (MGEX) and International GNSS Monitoring and Assessment System (iGMAS) are used in positioning experiments. The distribution of MGEX and iGMAS stations relating to the BDS-3 observations is presented in Figure2.Figure 2
Distribution of BDS-3 satellite tracking stations.


### 3.1. Analysis of BDS-3 Code Multipath Delays

Four groups of experiments are conducted with different schemes to analyze the proposed method. In the experiments, one-month observations (day of year (DOY) 214–244, 2020) made by BDS-3 stations are selected. The characteristics of the multipath delays of BDS-3 satellite code observations at different frequencies (i.e., B1I, B3I, B1C, and B2a) are first analyzed. Results for two iGMAS stations calculated using equations (1)–(7) are shown in Figure 3, where the series of multipath delays of station-WUH1 code observations of three satellites (i.e., C34, C38, and C59) and the corresponding satellite elevations are plotted for DOY 228 as examples. Figure 3 suggests that BDS-3 observations contain appreciable code multipath delays. Meanwhile, a trend term is seen along the satellite elevation, with the observations corresponding to lower elevations having more obvious multipath errors. Therefore, the trend term (systematic part) and the residual term (random part) should be modeled separately in the mitigation of BDS-3 multipath errors.Figure 3
Series of BDS-3 multipath delays for the iGMAS (WUH1) station (DOY 228).
(a)
C34 (MEO)(b)
C38 (IGSO)(c)
C59 (GEO)The first group of experiments are conducted to test the algorithm of sparse modeling in processing the BDS-3 code multipath series. In the experiments, three comparison schemes are designed to fully analyze the performance of the sparse modeling of the multipath delays of BDS-3 code observations.Scheme 1: A series of estimated BDS-3 multipath delays at different frequencies is first modeled using a regularization algorithm [20]. Typical distributions for MP1 (station MP at the B1I frequency) after regularization are plotted in Figure 4(b), where all medium Earth orbit (MEO) satellite observations of WUH1 on DOY 228 are input in the analysis.Figure 4
Typical distributions of the BDS-3 MP1 series obtained using the sparse modeling algorithm.
(a)
Original series(b)
Scheme 1(c)
Scheme 2(d)
Scheme 3Scheme 2: Similarly, the series of multipath delays are optimized using EN of the sparse structure. The processed series are presented in Figure4(c).The objective function of Scheme 2 is further improved by introducing the TV term. The corresponding typical distributions of the MP1 series are shown in Figure4(d).Owing to the large volume of experimental observations, the MEO satellites of BDS-3 observations are taken as examples to show the performances of the sparse modeling algorithm in Figure4. Furthermore, the results for frequencies B3I, B1C, and B2a are analyzed using the one-month observations for all BDS-3 MEO satellites, where the standard deviations of the series are given in Table 1. The table shows that different values of the standard deviation are output for different schemes, where the adoption of regularization, EN, and TV denoises the MP series, realizes sparse modeling, and overcomes the restriction of limited series, respectively.Table 1
Standard deviations (m) of multipath series at frequencies B1I, B3I, B1C, and B2a for BDS-3 satellite observations.
B1IB3IB1CB2aMEOIGSOMEOIGSOMEOIGSOMEOIGSOMGEXOriginal series0.3000.2950.2090.2870.2280.2250.2400.194Scheme 10.2860.2810.1990.2730.2160.2140.2280.185Scheme 20.1940.1930.1050.1910.1260.1330.1380.102Scheme 30.1690.1720.0970.1880.1140.1190.1290.098iGMASOriginal series0.3110.4070.2660.3680.3160.4950.2590.347Scheme 10.2960.3880.2540.3500.3010.4720.2470.331Scheme 20.2140.3330.1860.3030.2080.3670.1750.226Scheme 30.2090.2850.1800.2420.1960.3550.1710.224The results of multipath delays obtained using the different processing strategies indicate that sparse modeling with the aid of the TV term provides a more stable series. Because the series of observations selected in constructing the mitigation model of multipath delays are usually within a period of 1 day, the TV term added to the sparse modeling function theoretically reduces the effects of noise and limited epochs. An integrated modeling strategy of the BDS-3 multipath series is investigated in the second group of experiments. Three modeling scenarios of multipath delays are considered for four iGMAS and four MGEX stations. The details of integrated modeling were presented in Figure1. The experiments on multipath modeling using the proposed methods are described as follows.


#### 3.1.1. Scenario 1

On the basis of the results of Scheme 3, the series of BDS-3 multipath delays is firstly modeled using the AR algorithm. Model residuals are output and analyzed.


#### 3.1.2. Scenario 2

In mitigating the trend terms of the residuals ofScenario 1, a polynomial function is used to reduce the multipath errors. Meanwhile, the strategy of Scenario 1 is adopted to model residuals after the processing of trend terms.


#### 3.1.3. Scenario 3

With the consideration of the two-step strategy adopted inScenario 2, the integrated estimation of model coefficients is tested using the integrated LS-combined AR method.The model residuals of the three scenarios are shown in Figure5 for the analysis of the results of the integrated method in modeling the BDS-3 multipath series. The figure shows that the integrated method, especially in one-step estimation, is more stable than the traditional AR-only method. Therefore, the integrated estimation strategy of modeling multipath delays using the LS-combined AR method is recommended for the processing of BDS-3 code observations.Figure 5
Model residuals of BDS-3 B1I and B3I code observations (original series and post modeling) forScenarios 1, 2, and 3 (C33, C34, and C35).
(a)
AR method only (B1C)(b)
AR method only (B2a)(c)
Two-step strategy (B1C)(d)
Two-step strategy (B2a)(e)
Integrated LS-combined AR method (B1C)(f)
Integrated LS-combined AR method (B2a)Figure5 shows that the LS-combined AR-integrated strategy well models and mitigates the BDS-3 code multipath errors and can be induced to improve the quality of BDS-3 satellite observations. However, the real-time applications of BDS-3 services cannot estimate model parameters with long-term observations. As an example, widely adopted sidereal filtering is conducted on the basis of the satellite orbit period. Therefore, in correcting for the multipath delays of next-arc observations in real-time or near real-time applications, it is necessary to further improve the model coefficients of multipath errors. According to Figure 1, the multipath delays are output and modeled by LS-combined AR method based on equations (1)–(15) for BDS-3 satellites’ multifrequency observations. As mentioned above, we take the estimated model coefficients as prior information into the solution of multipath errors of the next-arc observations for real-time or rapid applications. In experiments, the tracking station of iGMAS (WUH1, DOY 228, 2020) is taken into modelling and correcting BDS-3 (MEO satellites) multipath delays. It should be noted that the proposed method takes the correlation between multipath delays and satellites’ elevations into consideration, which is different from the traditional sidereal filtering based on the repeat time of satellites’ orbit. Therefore, only one-day (DOY 228, 2020) observations are selected into experiments to model multipath delay and real-time correction with its prior information. In verifying the availability of the improved multipath model with the prior constraints, two experiments are conducted to analyze the multipath residuals. Additionally, the observations are preprocessed to detect and repair the cycle slips and delete satellites without observations during modeling multipath delays.


#### 3.1.4. Experiment 1

On the basis of the estimated model coefficients of a multipath series, the next-arc observations are corrected using the predicted values. The residuals between the predicted and real multipath series are plotted in Figure6(b) (AR method only) and Figure 6(c) (integrated-LS combined AR method).Figure 6
Residuals between predicted and estimated delays ofExperiments 1 and 2 in the third group of experiments (WUH1, DOY 228, 2020).
(a)
Original multipath series(b)
Multipath residuals for the AR method only(c)
Multipath residuals of the LS-combined AR(d)
Multipath residuals with prior information


#### 3.1.5. Experiment 2

In the correction of next-arc BDS-3 multipath errors, the estimated model coefficients are set as prior information to obtain the new model. The next-arc residuals between the real and calculated values are shown in Figure6(d).The residuals of the multipath series in the third group of experiments show that the next-arc BDS-3 observations can be optimized with the prior information in the estimation of multipath models. Table2 gives the average daily root-mean-square (RMS) values of multipath residuals for different frequencies (MGEX stations) to show the performance of the proposed strategy in Experiments 1 and 2. Figure 6 and Table 2 show that the improved strategy of modeling BDS-3 multipath delays indeed improves the accuracy of BDS-3 observations for rapid and real-time applications. In previous researches [2], the BDS-2 satellites’ code biases are discussed and corrected by a proposed method called the one-step method, in which the high-correction factor between multipath delays and code biases are considered. In this research, the multipath delays of BDS-3 multifrequency observations are analyzed, in which a systematic term and a random term divided from BDS-3 satellites’ code observations can be found from the first group of experiments. It should be noted that the BDS-2 code biases are different from BDS-3 systematic multipath delays based on mechanism, characteristic, and order of magnitude. Although the code bias related to the BDS-2 internal multipath vanished for BDS-3 code observations, the multipath of BDS-3 multifrequency observations should be paid more attention and be improved to optimize the positioning services. Due to the differences between BDS-2 and BDS-3 satellites, we usually take additional parameters, such ISB and DCB, to accurately estimate BDS-2- and BDS-3-related results [3]. In this research, it is difficult to mitigate the multipath delays of BDS-2 and BDS-3 within one solution, which reasons can be summarized as two points: firstly, differences of the coefficient of systematic terms for the BDS-2 and BDS-3 multipath delays are more than an order of magnitude [16]. So, we cannot together estimate the coefficients with a unified model; secondly, the code bias of BDS-2 observations should be firstly corrected. And the residuals of multipath of BDS-2 can be modeled with a same strategy of BDS-3. In the next step, we will take the proposed model into the positioning services of combined BDS-2 and BDS-3 satellites after the model and elimination of BDS-2 code biases.Table 2
Average daily RMS values (m) of multipath residuals in the third group of experiments.
B1IB3IB1CB2aExperiment 1 (LS combined AR)0.2610.1440.1720.154Experiment 2 (prior information)0.1350.1020.0920.089


### 3.2. Results for BDS-3 PPP Solutions

The aim of this study is at improving the performance of BDS-3 PPP services. It is concluded from the above three groups of experiments that the proposed integrated method with the aid of sparse modeling provides a more stable series of the BDS-3 multipath delay. Meanwhile, a fourth group of experiments on BDS-3 single-frequency PPP and double-frequency real-time PPP is conducted using the proposed strategy for the BDS-3 multipath delay. Further details of configurations for single-frequency PPP and real-time PPP have been published [2]. We set the cutoff angle at 5° in estimating the positioning parameters. The antenna types are described in the literature [3]. Because of the lack of accurate BDS-3 satellite phase center variation (PCV) parameters, we simply ignore the effects of PCV on the PPP solution. All positioning experiments are conducted using the orbit determination and positioning software package of CUM (an iGMAS analysis center), where the reference frame is ITRF 2014.


#### 3.2.1. Single-Frequency PPP Experiments

Three iGMAS stations and one MGEX station making BDS-3 observations are selected in the experiments for the analysis of BDS-3 single-frequency PPP. The configuration of single-frequency PPP is listed in Table3. The double-frequency combination is adopted to model multipath delays. The experimental session is set as DOY 221–228, 2020, for which observations made at four frequencies (i.e., BI1, B3I, B1C, and B2a) are analyzed. Precise satellite products, such as orbit and clock offsets, are downloaded from the GFZ analysis center (GBM products). Moreover, four schemes using different strategies to mitigate BDS-3 multipath delays, namely, the original series (none), AR model (AR), LS+AR (two-step), and integrated correction (one step), are analyzed. Residuals of station coordinates are plotted in Figure 7 (B1I) and Figure 8 (B1C), where the results for DOY 226–228 are presented with 2880 epochs per day. Moreover, to present all results of the four stations and four frequencies, the average RMS values of eight-day single-frequency PPP and the corresponding improvements are given in Table 4 (B1I and B3I) and Table 5 (B1C and B2a). In particular, the processing of different strategies for the single-frequency PPP solution is analyzed after the sparse modeling of BDS-3 multipath delays with the aid of the EN structure and TV term.Table 3
Strategy for the single-frequency PPP solution [2].
ParametersMain configurationsObservationsFour frequencies of phase and pseudorange (B1I, B3I, B1C, and B2a)Cutoff elevation angle5°Weigh assignmentElevation-dependent (below 30° by 1/2sin (elevation))Epoch interval30 seconds per epochIonospheric parametersPriorproducts+random walkTropospheric parametersSaastamonient+estimated wet componentsSatellites PCOhttp://www.igmas.orgTide modelsInternational earth rotation and reference systems (IERS) 2010RelativityIERS 2010Gravity modelEGM 08 model with12×12 ordersSatellites PCVIgnoredStation positionsEstimated and outputReceivers clock errorsEstimated and outputAmbiguity parametersFloat solutionsFigure 7
Results of BDS-3 single-frequency PPP obtained for B1I observations (DOY 226-228).Figure 8
Results of BDS-3 single-frequency PPP for B1C observations (DOY 226-228).Table 4
Average daily RMS values (m) of position residuals (B1I, B3I) for four stations obtained using three BDS-3 observation multipath delay processing strategies and the corresponding improvements.
StationsDirectionsB1IB3INoneTwo-stepImp. (%)One-stepImp. (%)NoneTwo-stepImp. (%)One-stepImp. (%)WUH2E−0.085−0.03756.5−0.03262.4−0.253−0.2491.6−0.13247.8N−0.044−0.0416.8−0.03325.0−0.091−0.0847.6−0.05638.5U−0.089−0.0845.6−0.05933.7−0.239−0.20514.2−0.15037.2KNDYE−0.024−0.01920.8−0.01825.0−0.042−0.03711.9−0.03419.0N−0.022−0.01340.9−0.01245.4−0.042-0.03614.3−0.02833.3U0.2590.22513.10.22214.30.1780.1751.60.1713.9BJF1E−0.049−0.0482.0−0.0474.1−0.045−0.0450.0−0.0426.6N−0.031−0.0296.5−0.0289.6−0.038−0.0365.2−0.0365.3U0.0580.05112.10.04325.90.033−0.0330.00.0323.0XIA1E−0.043−0.0399.3−0.03616.3−0.041−0.0394.9−0.0387.3N−0.029−0.02610.3−0.02417.2−0.045−0.0442.2−0.0426.6U0.5800.5770.50.5731.20.8860.8850.10.8840.2Table 5
Average daily RMS values (m) of position residuals (B1C, B2a) for four stations obtained using three BDS-3 observation multipath delay processing strategies and the corresponding improvements.
StationsDirectionsB1CB2aNoneTwo-stepImp. (%)One-stepImp. (%)NoneTwo-stepImp. (%)One-stepImp. (%)WUH2E−0.058−0.05013.8−0.04522.4−0.029−0.02513.8−0.02513.8N−0.038−0.03410.50.03215.80.0820.05236.60.05039.0U−0.249−0.2288.4−0.22210.8−0.271−0.24011.4−0.22218.1KNDYE−0.0500.01864.00.00982.0−0.115−0.08030.4−0.03569.6N−0.030−0.0293.3−0.01743.3−0.054−0.0523.7−0.0507.4U−0.105−0.1050.0−0.1013.8−0.240−0.07668.3−0.05975.4BJF1E−0.022−0.0220.0−0.01340.9−0.055−0.02849.1−0.02750.9N−0.028−0.0280.0−0.02414.3−0.017−0.015011.8−0.05070.6U−0.320−0.3141.9−0.28211.9−0.608−0.5922.6−0.07387.9XIA1E−0.032−0.0299.4−0.02812.5−0.017−0.0165.9−0.0165.9N−0.020−0.01620.0−0.01525.0−0.024−0.0234.2−0.0228.3U−0.075−0.0715.3−0.0698.0−0.033−0.02136.4−0.01457.6In the experiments, to avoid errors common to the different PPP solutions, positions of the double-frequency BDS/global positioning system (GPS) solution (i.e., the millimeter difference with the SINEX file) are taken as references to calculate the accuracy of single-frequency PPP. Figures7 and 8 show that the quality is obviously better for B1C than for B1I, which is consistent with the results of previous research [16]. Additionally, the integrated LS-combined AR strategy obtains position solutions more accurately. Specifically, the results of single-frequency PPP obtained for the four frequencies reveal that the accuracy in eastward (E), northward (N), and upward (U) directions is improved by at least 4.1%, 9.6%, and 1.2%, respectively, for B1I; 6.6%, 5.3%, and 0.2%, respectively, for B3I, 12.5%, 14.3%, and 3.8%, respectively, for B1C; and 5.9%, 7.4%, and 18.1%, respectively, for B2a, relative to the results obtained using the traditional method (none). However, because of the high-quality observation of BDS-3, the integrated strategy is slightly more optimized than the two-step method. The ultimate goal of the present study is the development of rapid and high-accuracy positioning services of BDS-3. Therefore, real-time double-frequency BDS-3 PPP is tested with different strategies for multipath errors in the following experiments.The real-time PPP solution is obtained according to the stage of CUM. A flowchart of real-time BDS-3 double-frequency PPP has been published in the literature [2]. Before obtaining the solution of real-time PPP, the satellite ultrarapid orbit and clock offsets are first estimated [2], which is not the focus of the present study. Additionally, the experimental session of real-time PPP is set as DOY 221–228, 2020, where the combinations of B1I+B3I and B1C+B2a are tested for four iGMAS and four MGEX stations. To avoid the restrictions of the real-time stream in the experiments, the RINEX files of observations are input into real-time PPP solutions. Similarly, four schemes of comparative experiments on BDS-3 multipath delays are designed, namely, no correction (none), correction with the AR model (AR), correction with the two-step LS+AR model (LS+AR), and the integrated estimation of the LS+AR model (one step). However, in comparison with BDS-3 single-frequency PPP, the AR strategy is adopted for the estimated coefficients of the previous day and the multipath delays are mitigated by the predicted values. In contrast, the one-step and two-step strategies are based on the prior information of previous model coefficients. Furthermore, the updated model coefficients are estimated with short-term observations (of an hour or so).The position residuals (BJF1) for three continuous days (DOY 226–228, 2020) are plotted in Figure9 (B1I+B3I) and Figure 10 (B1C+B2a) to show the results of BDS-3 real-time PPP. Moreover, the average daily RMS residuals for eight days, four stations, and three strategies (none, two-step, and one-step) and the corresponding improvements are listed in Table 6. Figures 9 and 10 show that there are no appreciable differences in the final precision among the four strategies, although the one-step correction performs slightly better that the two-step strategy. Furthermore, there is a bias in the direction of U, which might be due to the uncalibrated PCV and phase center offset of the satellites.Figure 9
Position residuals of BDS-3 double-frequency real-time PPP forB1I+B3I (BJF1).Figure 10
Position residuals of BDS-3 double-frequency real-time PPP forB1C+B2a (BJF1).Table 6
Average daily RMS values (m) of position residuals for four stations and three multipath error modeling strategies and corresponding improvements.
StationsDirectionsB1I+B3IB1C+B2aNoneTwo-stepImp. (%)One-stepImp. (%)NoneTwo-stepImp. (%)One-stepImp. (%)BJF1E−0.054−0.04025.9−0.03927.8−0.056−0.04814.3−0.04126.8N−0.040−0.02342.5−0.01660.0−0.043−0.03225.6−0.02346.5U−0.174−0.1579.8−0.09445.9−0.156−0.12718.6−0.09241.0KNDYE−0.055−0.0550.0−0.04321.8−0.047−0.0438.5−0.03819.1N−0.040−0.01660.0−0.01367.5−0.040−0.02537.5−0.02147.5U−0.204−0.09354.4−0.04378.9−0.161−0.1497.5−0.10336.0WUH1E−0.043−0.0430.0−0.03713.9−0.039−0.02633.3−0.02633.3N−0.035−0.01654.3−0.01362.9−0.039−0.02243.6−0.01756.4U−0.183−0.10045.4−0.04674.9−0.154−0.13114.9−0.09240.3XIA1E−0.052−0.01473.1−0.01375.0−0.051−0.0501.9−0.03433.3N−0.054−0.01572.2−0.01375.9−0.033−0.02815.2−0.01651.5U−0.108−0.05846.3−0.02874.1−0.142−0.1345.6−0.10823.9There are no obvious differences in the simulated final position residuals between single-frequency PPP and real-time PPP. Table6 shows that the improvements in E, N, and U directions are at least 13.9%, 60.0%, and 45.9%, respectively, and 19.1%, 46.5%, and 23.9%, respectively, for B1I+B3I and B1C+B2a, respectively. It is suggested that the one-step strategy with the prior information for modeling the BDS-3 multipath has better real-time application than the traditional method (none). It should be noted that the convergence time is also an important indicator for the PPP performances. In experiments, the convergence time and position residuals are simultaneously considered. However, it is suggested that there is an average improvement within 2% for E and N directions based on the proposed method, which cannot be as a significant improvement in experiments. However, it is difficult to count the convergence time of the U direction as a systematic error found. Therefore, the convergence time of the PPP solution is ignored in experiments.More BDS-3 observations collected by GNSS receivers are used to test the proposed strategy of mitigating the BDS-3 multipath delay. Furthermore, the results of the SPP experiment are verified using a set of real observations, the trajectory of which is shown in Figure11. In consideration of the decoding ability of GNSS receivers, the code observations of B1C are analyzed with an interval of 0.2 s, where B1I and B3I are used to model the code multipath delays. In experiments, three test scenarios are conducted with the same observations, strategies, and parameter estimations. However, the BDS-3 code multipath delays are modeled and corrected with the proposed method (one-step), traditional method (two-step), and none considerations (none) in the preprocessing of BDS-3 observations. In a word, the changes of position residuals are caused by the different multipath mitigation method. The single-epoch solution of SPP in E and N directions is shown for three strategies (none, one-step, and two-step) in Figure 12. In analyzing the accuracy of SPP, the station positions of GPS double-frequency kinematic PPP are taken as a reference. SPP does not provide the desired result owing to the complex kinematic environments around the stations. However, the precision and stability can be improved easily by adopting the one-step and two-step strategies. Moreover, there are no appreciable differences in results between the one-step and two-step strategies as no prior information is used in modeling BDS-3 multipath errors. Therefore, a more precise and accurate strategy should be developed to outperform the proposed one-step correction for multipath delays in real-time applications. Figure 6 shows that the multiple paths resulting from the kinematic environment cannot be precisely corrected for, whereas the systematic and trend terms along the elevation angle can be reduced using our method. Because of the complex environments, the improvements of dynamic multipath mitigation for position residuals cannot totally be ascribed to the proposed method as the significant errors in different directions. In the next step, the dynamic multipath mitigation will be paid more attention to improve the applications of BDS-3 real-time services.Figure 11
Trajectory of the SPP experiment.Figure 12
Position residuals obtained in SPP experiments with different strategies.


## 3.1. Analysis of BDS-3 Code Multipath Delays

Four groups of experiments are conducted with different schemes to analyze the proposed method. In the experiments, one-month observations (day of year (DOY) 214–244, 2020) made by BDS-3 stations are selected. The characteristics of the multipath delays of BDS-3 satellite code observations at different frequencies (i.e., B1I, B3I, B1C, and B2a) are first analyzed. Results for two iGMAS stations calculated using equations (1)–(7) are shown in Figure 3, where the series of multipath delays of station-WUH1 code observations of three satellites (i.e., C34, C38, and C59) and the corresponding satellite elevations are plotted for DOY 228 as examples. Figure 3 suggests that BDS-3 observations contain appreciable code multipath delays. Meanwhile, a trend term is seen along the satellite elevation, with the observations corresponding to lower elevations having more obvious multipath errors. Therefore, the trend term (systematic part) and the residual term (random part) should be modeled separately in the mitigation of BDS-3 multipath errors.Figure 3
Series of BDS-3 multipath delays for the iGMAS (WUH1) station (DOY 228).
(a)
C34 (MEO)(b)
C38 (IGSO)(c)
C59 (GEO)The first group of experiments are conducted to test the algorithm of sparse modeling in processing the BDS-3 code multipath series. In the experiments, three comparison schemes are designed to fully analyze the performance of the sparse modeling of the multipath delays of BDS-3 code observations.Scheme 1: A series of estimated BDS-3 multipath delays at different frequencies is first modeled using a regularization algorithm [20]. Typical distributions for MP1 (station MP at the B1I frequency) after regularization are plotted in Figure 4(b), where all medium Earth orbit (MEO) satellite observations of WUH1 on DOY 228 are input in the analysis.Figure 4
Typical distributions of the BDS-3 MP1 series obtained using the sparse modeling algorithm.
(a)
Original series(b)
Scheme 1(c)
Scheme 2(d)
Scheme 3Scheme 2: Similarly, the series of multipath delays are optimized using EN of the sparse structure. The processed series are presented in Figure4(c).The objective function of Scheme 2 is further improved by introducing the TV term. The corresponding typical distributions of the MP1 series are shown in Figure4(d).Owing to the large volume of experimental observations, the MEO satellites of BDS-3 observations are taken as examples to show the performances of the sparse modeling algorithm in Figure4. Furthermore, the results for frequencies B3I, B1C, and B2a are analyzed using the one-month observations for all BDS-3 MEO satellites, where the standard deviations of the series are given in Table 1. The table shows that different values of the standard deviation are output for different schemes, where the adoption of regularization, EN, and TV denoises the MP series, realizes sparse modeling, and overcomes the restriction of limited series, respectively.Table 1
Standard deviations (m) of multipath series at frequencies B1I, B3I, B1C, and B2a for BDS-3 satellite observations.
B1IB3IB1CB2aMEOIGSOMEOIGSOMEOIGSOMEOIGSOMGEXOriginal series0.3000.2950.2090.2870.2280.2250.2400.194Scheme 10.2860.2810.1990.2730.2160.2140.2280.185Scheme 20.1940.1930.1050.1910.1260.1330.1380.102Scheme 30.1690.1720.0970.1880.1140.1190.1290.098iGMASOriginal series0.3110.4070.2660.3680.3160.4950.2590.347Scheme 10.2960.3880.2540.3500.3010.4720.2470.331Scheme 20.2140.3330.1860.3030.2080.3670.1750.226Scheme 30.2090.2850.1800.2420.1960.3550.1710.224The results of multipath delays obtained using the different processing strategies indicate that sparse modeling with the aid of the TV term provides a more stable series. Because the series of observations selected in constructing the mitigation model of multipath delays are usually within a period of 1 day, the TV term added to the sparse modeling function theoretically reduces the effects of noise and limited epochs. An integrated modeling strategy of the BDS-3 multipath series is investigated in the second group of experiments. Three modeling scenarios of multipath delays are considered for four iGMAS and four MGEX stations. The details of integrated modeling were presented in Figure1. The experiments on multipath modeling using the proposed methods are described as follows.


### 3.1.1. Scenario 1

On the basis of the results of Scheme 3, the series of BDS-3 multipath delays is firstly modeled using the AR algorithm. Model residuals are output and analyzed.


### 3.1.2. Scenario 2

In mitigating the trend terms of the residuals ofScenario 1, a polynomial function is used to reduce the multipath errors. Meanwhile, the strategy of Scenario 1 is adopted to model residuals after the processing of trend terms.


### 3.1.3. Scenario 3

With the consideration of the two-step strategy adopted inScenario 2, the integrated estimation of model coefficients is tested using the integrated LS-combined AR method.The model residuals of the three scenarios are shown in Figure5 for the analysis of the results of the integrated method in modeling the BDS-3 multipath series. The figure shows that the integrated method, especially in one-step estimation, is more stable than the traditional AR-only method. Therefore, the integrated estimation strategy of modeling multipath delays using the LS-combined AR method is recommended for the processing of BDS-3 code observations.Figure 5
Model residuals of BDS-3 B1I and B3I code observations (original series and post modeling) forScenarios 1, 2, and 3 (C33, C34, and C35).
(a)
AR method only (B1C)(b)
AR method only (B2a)(c)
Two-step strategy (B1C)(d)
Two-step strategy (B2a)(e)
Integrated LS-combined AR method (B1C)(f)
Integrated LS-combined AR method (B2a)Figure5 shows that the LS-combined AR-integrated strategy well models and mitigates the BDS-3 code multipath errors and can be induced to improve the quality of BDS-3 satellite observations. However, the real-time applications of BDS-3 services cannot estimate model parameters with long-term observations. As an example, widely adopted sidereal filtering is conducted on the basis of the satellite orbit period. Therefore, in correcting for the multipath delays of next-arc observations in real-time or near real-time applications, it is necessary to further improve the model coefficients of multipath errors. According to Figure 1, the multipath delays are output and modeled by LS-combined AR method based on equations (1)–(15) for BDS-3 satellites’ multifrequency observations. As mentioned above, we take the estimated model coefficients as prior information into the solution of multipath errors of the next-arc observations for real-time or rapid applications. In experiments, the tracking station of iGMAS (WUH1, DOY 228, 2020) is taken into modelling and correcting BDS-3 (MEO satellites) multipath delays. It should be noted that the proposed method takes the correlation between multipath delays and satellites’ elevations into consideration, which is different from the traditional sidereal filtering based on the repeat time of satellites’ orbit. Therefore, only one-day (DOY 228, 2020) observations are selected into experiments to model multipath delay and real-time correction with its prior information. In verifying the availability of the improved multipath model with the prior constraints, two experiments are conducted to analyze the multipath residuals. Additionally, the observations are preprocessed to detect and repair the cycle slips and delete satellites without observations during modeling multipath delays.


### 3.1.4. Experiment 1

On the basis of the estimated model coefficients of a multipath series, the next-arc observations are corrected using the predicted values. The residuals between the predicted and real multipath series are plotted in Figure6(b) (AR method only) and Figure 6(c) (integrated-LS combined AR method).Figure 6
Residuals between predicted and estimated delays ofExperiments 1 and 2 in the third group of experiments (WUH1, DOY 228, 2020).
(a)
Original multipath series(b)
Multipath residuals for the AR method only(c)
Multipath residuals of the LS-combined AR(d)
Multipath residuals with prior information


### 3.1.5. Experiment 2

In the correction of next-arc BDS-3 multipath errors, the estimated model coefficients are set as prior information to obtain the new model. The next-arc residuals between the real and calculated values are shown in Figure6(d).The residuals of the multipath series in the third group of experiments show that the next-arc BDS-3 observations can be optimized with the prior information in the estimation of multipath models. Table2 gives the average daily root-mean-square (RMS) values of multipath residuals for different frequencies (MGEX stations) to show the performance of the proposed strategy in Experiments 1 and 2. Figure 6 and Table 2 show that the improved strategy of modeling BDS-3 multipath delays indeed improves the accuracy of BDS-3 observations for rapid and real-time applications. In previous researches [2], the BDS-2 satellites’ code biases are discussed and corrected by a proposed method called the one-step method, in which the high-correction factor between multipath delays and code biases are considered. In this research, the multipath delays of BDS-3 multifrequency observations are analyzed, in which a systematic term and a random term divided from BDS-3 satellites’ code observations can be found from the first group of experiments. It should be noted that the BDS-2 code biases are different from BDS-3 systematic multipath delays based on mechanism, characteristic, and order of magnitude. Although the code bias related to the BDS-2 internal multipath vanished for BDS-3 code observations, the multipath of BDS-3 multifrequency observations should be paid more attention and be improved to optimize the positioning services. Due to the differences between BDS-2 and BDS-3 satellites, we usually take additional parameters, such ISB and DCB, to accurately estimate BDS-2- and BDS-3-related results [3]. In this research, it is difficult to mitigate the multipath delays of BDS-2 and BDS-3 within one solution, which reasons can be summarized as two points: firstly, differences of the coefficient of systematic terms for the BDS-2 and BDS-3 multipath delays are more than an order of magnitude [16]. So, we cannot together estimate the coefficients with a unified model; secondly, the code bias of BDS-2 observations should be firstly corrected. And the residuals of multipath of BDS-2 can be modeled with a same strategy of BDS-3. In the next step, we will take the proposed model into the positioning services of combined BDS-2 and BDS-3 satellites after the model and elimination of BDS-2 code biases.Table 2
Average daily RMS values (m) of multipath residuals in the third group of experiments.
B1IB3IB1CB2aExperiment 1 (LS combined AR)0.2610.1440.1720.154Experiment 2 (prior information)0.1350.1020.0920.089


## 3.1.1. Scenario 1

On the basis of the results of Scheme 3, the series of BDS-3 multipath delays is firstly modeled using the AR algorithm. Model residuals are output and analyzed.


## 3.1.2. Scenario 2

In mitigating the trend terms of the residuals ofScenario 1, a polynomial function is used to reduce the multipath errors. Meanwhile, the strategy of Scenario 1 is adopted to model residuals after the processing of trend terms.


## 3.1.3. Scenario 3

With the consideration of the two-step strategy adopted inScenario 2, the integrated estimation of model coefficients is tested using the integrated LS-combined AR method.The model residuals of the three scenarios are shown in Figure5 for the analysis of the results of the integrated method in modeling the BDS-3 multipath series. The figure shows that the integrated method, especially in one-step estimation, is more stable than the traditional AR-only method. Therefore, the integrated estimation strategy of modeling multipath delays using the LS-combined AR method is recommended for the processing of BDS-3 code observations.Figure 5
Model residuals of BDS-3 B1I and B3I code observations (original series and post modeling) forScenarios 1, 2, and 3 (C33, C34, and C35).
(a)
AR method only (B1C)(b)
AR method only (B2a)(c)
Two-step strategy (B1C)(d)
Two-step strategy (B2a)(e)
Integrated LS-combined AR method (B1C)(f)
Integrated LS-combined AR method (B2a)Figure5 shows that the LS-combined AR-integrated strategy well models and mitigates the BDS-3 code multipath errors and can be induced to improve the quality of BDS-3 satellite observations. However, the real-time applications of BDS-3 services cannot estimate model parameters with long-term observations. As an example, widely adopted sidereal filtering is conducted on the basis of the satellite orbit period. Therefore, in correcting for the multipath delays of next-arc observations in real-time or near real-time applications, it is necessary to further improve the model coefficients of multipath errors. According to Figure 1, the multipath delays are output and modeled by LS-combined AR method based on equations (1)–(15) for BDS-3 satellites’ multifrequency observations. As mentioned above, we take the estimated model coefficients as prior information into the solution of multipath errors of the next-arc observations for real-time or rapid applications. In experiments, the tracking station of iGMAS (WUH1, DOY 228, 2020) is taken into modelling and correcting BDS-3 (MEO satellites) multipath delays. It should be noted that the proposed method takes the correlation between multipath delays and satellites’ elevations into consideration, which is different from the traditional sidereal filtering based on the repeat time of satellites’ orbit. Therefore, only one-day (DOY 228, 2020) observations are selected into experiments to model multipath delay and real-time correction with its prior information. In verifying the availability of the improved multipath model with the prior constraints, two experiments are conducted to analyze the multipath residuals. Additionally, the observations are preprocessed to detect and repair the cycle slips and delete satellites without observations during modeling multipath delays.


## 3.1.4. Experiment 1

On the basis of the estimated model coefficients of a multipath series, the next-arc observations are corrected using the predicted values. The residuals between the predicted and real multipath series are plotted in Figure6(b) (AR method only) and Figure 6(c) (integrated-LS combined AR method).Figure 6
Residuals between predicted and estimated delays ofExperiments 1 and 2 in the third group of experiments (WUH1, DOY 228, 2020).
(a)
Original multipath series(b)
Multipath residuals for the AR method only(c)
Multipath residuals of the LS-combined AR(d)
Multipath residuals with prior information


## 3.1.5. Experiment 2

In the correction of next-arc BDS-3 multipath errors, the estimated model coefficients are set as prior information to obtain the new model. The next-arc residuals between the real and calculated values are shown in Figure6(d).The residuals of the multipath series in the third group of experiments show that the next-arc BDS-3 observations can be optimized with the prior information in the estimation of multipath models. Table2 gives the average daily root-mean-square (RMS) values of multipath residuals for different frequencies (MGEX stations) to show the performance of the proposed strategy in Experiments 1 and 2. Figure 6 and Table 2 show that the improved strategy of modeling BDS-3 multipath delays indeed improves the accuracy of BDS-3 observations for rapid and real-time applications. In previous researches [2], the BDS-2 satellites’ code biases are discussed and corrected by a proposed method called the one-step method, in which the high-correction factor between multipath delays and code biases are considered. In this research, the multipath delays of BDS-3 multifrequency observations are analyzed, in which a systematic term and a random term divided from BDS-3 satellites’ code observations can be found from the first group of experiments. It should be noted that the BDS-2 code biases are different from BDS-3 systematic multipath delays based on mechanism, characteristic, and order of magnitude. Although the code bias related to the BDS-2 internal multipath vanished for BDS-3 code observations, the multipath of BDS-3 multifrequency observations should be paid more attention and be improved to optimize the positioning services. Due to the differences between BDS-2 and BDS-3 satellites, we usually take additional parameters, such ISB and DCB, to accurately estimate BDS-2- and BDS-3-related results [3]. In this research, it is difficult to mitigate the multipath delays of BDS-2 and BDS-3 within one solution, which reasons can be summarized as two points: firstly, differences of the coefficient of systematic terms for the BDS-2 and BDS-3 multipath delays are more than an order of magnitude [16]. So, we cannot together estimate the coefficients with a unified model; secondly, the code bias of BDS-2 observations should be firstly corrected. And the residuals of multipath of BDS-2 can be modeled with a same strategy of BDS-3. In the next step, we will take the proposed model into the positioning services of combined BDS-2 and BDS-3 satellites after the model and elimination of BDS-2 code biases.Table 2
Average daily RMS values (m) of multipath residuals in the third group of experiments.
B1IB3IB1CB2aExperiment 1 (LS combined AR)0.2610.1440.1720.154Experiment 2 (prior information)0.1350.1020.0920.089


## 3.2. Results for BDS-3 PPP Solutions

The aim of this study is at improving the performance of BDS-3 PPP services. It is concluded from the above three groups of experiments that the proposed integrated method with the aid of sparse modeling provides a more stable series of the BDS-3 multipath delay. Meanwhile, a fourth group of experiments on BDS-3 single-frequency PPP and double-frequency real-time PPP is conducted using the proposed strategy for the BDS-3 multipath delay. Further details of configurations for single-frequency PPP and real-time PPP have been published [2]. We set the cutoff angle at 5° in estimating the positioning parameters. The antenna types are described in the literature [3]. Because of the lack of accurate BDS-3 satellite phase center variation (PCV) parameters, we simply ignore the effects of PCV on the PPP solution. All positioning experiments are conducted using the orbit determination and positioning software package of CUM (an iGMAS analysis center), where the reference frame is ITRF 2014.


### 3.2.1. Single-Frequency PPP Experiments

Three iGMAS stations and one MGEX station making BDS-3 observations are selected in the experiments for the analysis of BDS-3 single-frequency PPP. The configuration of single-frequency PPP is listed in Table3. The double-frequency combination is adopted to model multipath delays. The experimental session is set as DOY 221–228, 2020, for which observations made at four frequencies (i.e., BI1, B3I, B1C, and B2a) are analyzed. Precise satellite products, such as orbit and clock offsets, are downloaded from the GFZ analysis center (GBM products). Moreover, four schemes using different strategies to mitigate BDS-3 multipath delays, namely, the original series (none), AR model (AR), LS+AR (two-step), and integrated correction (one step), are analyzed. Residuals of station coordinates are plotted in Figure 7 (B1I) and Figure 8 (B1C), where the results for DOY 226–228 are presented with 2880 epochs per day. Moreover, to present all results of the four stations and four frequencies, the average RMS values of eight-day single-frequency PPP and the corresponding improvements are given in Table 4 (B1I and B3I) and Table 5 (B1C and B2a). In particular, the processing of different strategies for the single-frequency PPP solution is analyzed after the sparse modeling of BDS-3 multipath delays with the aid of the EN structure and TV term.Table 3
Strategy for the single-frequency PPP solution [2].
ParametersMain configurationsObservationsFour frequencies of phase and pseudorange (B1I, B3I, B1C, and B2a)Cutoff elevation angle5°Weigh assignmentElevation-dependent (below 30° by 1/2sin (elevation))Epoch interval30 seconds per epochIonospheric parametersPriorproducts+random walkTropospheric parametersSaastamonient+estimated wet componentsSatellites PCOhttp://www.igmas.orgTide modelsInternational earth rotation and reference systems (IERS) 2010RelativityIERS 2010Gravity modelEGM 08 model with12×12 ordersSatellites PCVIgnoredStation positionsEstimated and outputReceivers clock errorsEstimated and outputAmbiguity parametersFloat solutionsFigure 7
Results of BDS-3 single-frequency PPP obtained for B1I observations (DOY 226-228).Figure 8
Results of BDS-3 single-frequency PPP for B1C observations (DOY 226-228).Table 4
Average daily RMS values (m) of position residuals (B1I, B3I) for four stations obtained using three BDS-3 observation multipath delay processing strategies and the corresponding improvements.
StationsDirectionsB1IB3INoneTwo-stepImp. (%)One-stepImp. (%)NoneTwo-stepImp. (%)One-stepImp. (%)WUH2E−0.085−0.03756.5−0.03262.4−0.253−0.2491.6−0.13247.8N−0.044−0.0416.8−0.03325.0−0.091−0.0847.6−0.05638.5U−0.089−0.0845.6−0.05933.7−0.239−0.20514.2−0.15037.2KNDYE−0.024−0.01920.8−0.01825.0−0.042−0.03711.9−0.03419.0N−0.022−0.01340.9−0.01245.4−0.042-0.03614.3−0.02833.3U0.2590.22513.10.22214.30.1780.1751.60.1713.9BJF1E−0.049−0.0482.0−0.0474.1−0.045−0.0450.0−0.0426.6N−0.031−0.0296.5−0.0289.6−0.038−0.0365.2−0.0365.3U0.0580.05112.10.04325.90.033−0.0330.00.0323.0XIA1E−0.043−0.0399.3−0.03616.3−0.041−0.0394.9−0.0387.3N−0.029−0.02610.3−0.02417.2−0.045−0.0442.2−0.0426.6U0.5800.5770.50.5731.20.8860.8850.10.8840.2Table 5
Average daily RMS values (m) of position residuals (B1C, B2a) for four stations obtained using three BDS-3 observation multipath delay processing strategies and the corresponding improvements.
StationsDirectionsB1CB2aNoneTwo-stepImp. (%)One-stepImp. (%)NoneTwo-stepImp. (%)One-stepImp. (%)WUH2E−0.058−0.05013.8−0.04522.4−0.029−0.02513.8−0.02513.8N−0.038−0.03410.50.03215.80.0820.05236.60.05039.0U−0.249−0.2288.4−0.22210.8−0.271−0.24011.4−0.22218.1KNDYE−0.0500.01864.00.00982.0−0.115−0.08030.4−0.03569.6N−0.030−0.0293.3−0.01743.3−0.054−0.0523.7−0.0507.4U−0.105−0.1050.0−0.1013.8−0.240−0.07668.3−0.05975.4BJF1E−0.022−0.0220.0−0.01340.9−0.055−0.02849.1−0.02750.9N−0.028−0.0280.0−0.02414.3−0.017−0.015011.8−0.05070.6U−0.320−0.3141.9−0.28211.9−0.608−0.5922.6−0.07387.9XIA1E−0.032−0.0299.4−0.02812.5−0.017−0.0165.9−0.0165.9N−0.020−0.01620.0−0.01525.0−0.024−0.0234.2−0.0228.3U−0.075−0.0715.3−0.0698.0−0.033−0.02136.4−0.01457.6In the experiments, to avoid errors common to the different PPP solutions, positions of the double-frequency BDS/global positioning system (GPS) solution (i.e., the millimeter difference with the SINEX file) are taken as references to calculate the accuracy of single-frequency PPP. Figures7 and 8 show that the quality is obviously better for B1C than for B1I, which is consistent with the results of previous research [16]. Additionally, the integrated LS-combined AR strategy obtains position solutions more accurately. Specifically, the results of single-frequency PPP obtained for the four frequencies reveal that the accuracy in eastward (E), northward (N), and upward (U) directions is improved by at least 4.1%, 9.6%, and 1.2%, respectively, for B1I; 6.6%, 5.3%, and 0.2%, respectively, for B3I, 12.5%, 14.3%, and 3.8%, respectively, for B1C; and 5.9%, 7.4%, and 18.1%, respectively, for B2a, relative to the results obtained using the traditional method (none). However, because of the high-quality observation of BDS-3, the integrated strategy is slightly more optimized than the two-step method. The ultimate goal of the present study is the development of rapid and high-accuracy positioning services of BDS-3. Therefore, real-time double-frequency BDS-3 PPP is tested with different strategies for multipath errors in the following experiments.The real-time PPP solution is obtained according to the stage of CUM. A flowchart of real-time BDS-3 double-frequency PPP has been published in the literature [2]. Before obtaining the solution of real-time PPP, the satellite ultrarapid orbit and clock offsets are first estimated [2], which is not the focus of the present study. Additionally, the experimental session of real-time PPP is set as DOY 221–228, 2020, where the combinations of B1I+B3I and B1C+B2a are tested for four iGMAS and four MGEX stations. To avoid the restrictions of the real-time stream in the experiments, the RINEX files of observations are input into real-time PPP solutions. Similarly, four schemes of comparative experiments on BDS-3 multipath delays are designed, namely, no correction (none), correction with the AR model (AR), correction with the two-step LS+AR model (LS+AR), and the integrated estimation of the LS+AR model (one step). However, in comparison with BDS-3 single-frequency PPP, the AR strategy is adopted for the estimated coefficients of the previous day and the multipath delays are mitigated by the predicted values. In contrast, the one-step and two-step strategies are based on the prior information of previous model coefficients. Furthermore, the updated model coefficients are estimated with short-term observations (of an hour or so).The position residuals (BJF1) for three continuous days (DOY 226–228, 2020) are plotted in Figure9 (B1I+B3I) and Figure 10 (B1C+B2a) to show the results of BDS-3 real-time PPP. Moreover, the average daily RMS residuals for eight days, four stations, and three strategies (none, two-step, and one-step) and the corresponding improvements are listed in Table 6. Figures 9 and 10 show that there are no appreciable differences in the final precision among the four strategies, although the one-step correction performs slightly better that the two-step strategy. Furthermore, there is a bias in the direction of U, which might be due to the uncalibrated PCV and phase center offset of the satellites.Figure 9
Position residuals of BDS-3 double-frequency real-time PPP forB1I+B3I (BJF1).Figure 10
Position residuals of BDS-3 double-frequency real-time PPP forB1C+B2a (BJF1).Table 6
Average daily RMS values (m) of position residuals for four stations and three multipath error modeling strategies and corresponding improvements.
StationsDirectionsB1I+B3IB1C+B2aNoneTwo-stepImp. (%)One-stepImp. (%)NoneTwo-stepImp. (%)One-stepImp. (%)BJF1E−0.054−0.04025.9−0.03927.8−0.056−0.04814.3−0.04126.8N−0.040−0.02342.5−0.01660.0−0.043−0.03225.6−0.02346.5U−0.174−0.1579.8−0.09445.9−0.156−0.12718.6−0.09241.0KNDYE−0.055−0.0550.0−0.04321.8−0.047−0.0438.5−0.03819.1N−0.040−0.01660.0−0.01367.5−0.040−0.02537.5−0.02147.5U−0.204−0.09354.4−0.04378.9−0.161−0.1497.5−0.10336.0WUH1E−0.043−0.0430.0−0.03713.9−0.039−0.02633.3−0.02633.3N−0.035−0.01654.3−0.01362.9−0.039−0.02243.6−0.01756.4U−0.183−0.10045.4−0.04674.9−0.154−0.13114.9−0.09240.3XIA1E−0.052−0.01473.1−0.01375.0−0.051−0.0501.9−0.03433.3N−0.054−0.01572.2−0.01375.9−0.033−0.02815.2−0.01651.5U−0.108−0.05846.3−0.02874.1−0.142−0.1345.6−0.10823.9There are no obvious differences in the simulated final position residuals between single-frequency PPP and real-time PPP. Table6 shows that the improvements in E, N, and U directions are at least 13.9%, 60.0%, and 45.9%, respectively, and 19.1%, 46.5%, and 23.9%, respectively, for B1I+B3I and B1C+B2a, respectively. It is suggested that the one-step strategy with the prior information for modeling the BDS-3 multipath has better real-time application than the traditional method (none). It should be noted that the convergence time is also an important indicator for the PPP performances. In experiments, the convergence time and position residuals are simultaneously considered. However, it is suggested that there is an average improvement within 2% for E and N directions based on the proposed method, which cannot be as a significant improvement in experiments. However, it is difficult to count the convergence time of the U direction as a systematic error found. Therefore, the convergence time of the PPP solution is ignored in experiments.More BDS-3 observations collected by GNSS receivers are used to test the proposed strategy of mitigating the BDS-3 multipath delay. Furthermore, the results of the SPP experiment are verified using a set of real observations, the trajectory of which is shown in Figure11. In consideration of the decoding ability of GNSS receivers, the code observations of B1C are analyzed with an interval of 0.2 s, where B1I and B3I are used to model the code multipath delays. In experiments, three test scenarios are conducted with the same observations, strategies, and parameter estimations. However, the BDS-3 code multipath delays are modeled and corrected with the proposed method (one-step), traditional method (two-step), and none considerations (none) in the preprocessing of BDS-3 observations. In a word, the changes of position residuals are caused by the different multipath mitigation method. The single-epoch solution of SPP in E and N directions is shown for three strategies (none, one-step, and two-step) in Figure 12. In analyzing the accuracy of SPP, the station positions of GPS double-frequency kinematic PPP are taken as a reference. SPP does not provide the desired result owing to the complex kinematic environments around the stations. However, the precision and stability can be improved easily by adopting the one-step and two-step strategies. Moreover, there are no appreciable differences in results between the one-step and two-step strategies as no prior information is used in modeling BDS-3 multipath errors. Therefore, a more precise and accurate strategy should be developed to outperform the proposed one-step correction for multipath delays in real-time applications. Figure 6 shows that the multiple paths resulting from the kinematic environment cannot be precisely corrected for, whereas the systematic and trend terms along the elevation angle can be reduced using our method. Because of the complex environments, the improvements of dynamic multipath mitigation for position residuals cannot totally be ascribed to the proposed method as the significant errors in different directions. In the next step, the dynamic multipath mitigation will be paid more attention to improve the applications of BDS-3 real-time services.Figure 11
Trajectory of the SPP experiment.Figure 12
Position residuals obtained in SPP experiments with different strategies.


## 3.2.1. Single-Frequency PPP Experiments

Three iGMAS stations and one MGEX station making BDS-3 observations are selected in the experiments for the analysis of BDS-3 single-frequency PPP. The configuration of single-frequency PPP is listed in Table3. The double-frequency combination is adopted to model multipath delays. The experimental session is set as DOY 221–228, 2020, for which observations made at four frequencies (i.e., BI1, B3I, B1C, and B2a) are analyzed. Precise satellite products, such as orbit and clock offsets, are downloaded from the GFZ analysis center (GBM products). Moreover, four schemes using different strategies to mitigate BDS-3 multipath delays, namely, the original series (none), AR model (AR), LS+AR (two-step), and integrated correction (one step), are analyzed. Residuals of station coordinates are plotted in Figure 7 (B1I) and Figure 8 (B1C), where the results for DOY 226–228 are presented with 2880 epochs per day. Moreover, to present all results of the four stations and four frequencies, the average RMS values of eight-day single-frequency PPP and the corresponding improvements are given in Table 4 (B1I and B3I) and Table 5 (B1C and B2a). In particular, the processing of different strategies for the single-frequency PPP solution is analyzed after the sparse modeling of BDS-3 multipath delays with the aid of the EN structure and TV term.Table 3
Strategy for the single-frequency PPP solution [2].
ParametersMain configurationsObservationsFour frequencies of phase and pseudorange (B1I, B3I, B1C, and B2a)Cutoff elevation angle5°Weigh assignmentElevation-dependent (below 30° by 1/2sin (elevation))Epoch interval30 seconds per epochIonospheric parametersPriorproducts+random walkTropospheric parametersSaastamonient+estimated wet componentsSatellites PCOhttp://www.igmas.orgTide modelsInternational earth rotation and reference systems (IERS) 2010RelativityIERS 2010Gravity modelEGM 08 model with12×12 ordersSatellites PCVIgnoredStation positionsEstimated and outputReceivers clock errorsEstimated and outputAmbiguity parametersFloat solutionsFigure 7
Results of BDS-3 single-frequency PPP obtained for B1I observations (DOY 226-228).Figure 8
Results of BDS-3 single-frequency PPP for B1C observations (DOY 226-228).Table 4
Average daily RMS values (m) of position residuals (B1I, B3I) for four stations obtained using three BDS-3 observation multipath delay processing strategies and the corresponding improvements.
StationsDirectionsB1IB3INoneTwo-stepImp. (%)One-stepImp. (%)NoneTwo-stepImp. (%)One-stepImp. (%)WUH2E−0.085−0.03756.5−0.03262.4−0.253−0.2491.6−0.13247.8N−0.044−0.0416.8−0.03325.0−0.091−0.0847.6−0.05638.5U−0.089−0.0845.6−0.05933.7−0.239−0.20514.2−0.15037.2KNDYE−0.024−0.01920.8−0.01825.0−0.042−0.03711.9−0.03419.0N−0.022−0.01340.9−0.01245.4−0.042-0.03614.3−0.02833.3U0.2590.22513.10.22214.30.1780.1751.60.1713.9BJF1E−0.049−0.0482.0−0.0474.1−0.045−0.0450.0−0.0426.6N−0.031−0.0296.5−0.0289.6−0.038−0.0365.2−0.0365.3U0.0580.05112.10.04325.90.033−0.0330.00.0323.0XIA1E−0.043−0.0399.3−0.03616.3−0.041−0.0394.9−0.0387.3N−0.029−0.02610.3−0.02417.2−0.045−0.0442.2−0.0426.6U0.5800.5770.50.5731.20.8860.8850.10.8840.2Table 5
Average daily RMS values (m) of position residuals (B1C, B2a) for four stations obtained using three BDS-3 observation multipath delay processing strategies and the corresponding improvements.
StationsDirectionsB1CB2aNoneTwo-stepImp. (%)One-stepImp. (%)NoneTwo-stepImp. (%)One-stepImp. (%)WUH2E−0.058−0.05013.8−0.04522.4−0.029−0.02513.8−0.02513.8N−0.038−0.03410.50.03215.80.0820.05236.60.05039.0U−0.249−0.2288.4−0.22210.8−0.271−0.24011.4−0.22218.1KNDYE−0.0500.01864.00.00982.0−0.115−0.08030.4−0.03569.6N−0.030−0.0293.3−0.01743.3−0.054−0.0523.7−0.0507.4U−0.105−0.1050.0−0.1013.8−0.240−0.07668.3−0.05975.4BJF1E−0.022−0.0220.0−0.01340.9−0.055−0.02849.1−0.02750.9N−0.028−0.0280.0−0.02414.3−0.017−0.015011.8−0.05070.6U−0.320−0.3141.9−0.28211.9−0.608−0.5922.6−0.07387.9XIA1E−0.032−0.0299.4−0.02812.5−0.017−0.0165.9−0.0165.9N−0.020−0.01620.0−0.01525.0−0.024−0.0234.2−0.0228.3U−0.075−0.0715.3−0.0698.0−0.033−0.02136.4−0.01457.6In the experiments, to avoid errors common to the different PPP solutions, positions of the double-frequency BDS/global positioning system (GPS) solution (i.e., the millimeter difference with the SINEX file) are taken as references to calculate the accuracy of single-frequency PPP. Figures7 and 8 show that the quality is obviously better for B1C than for B1I, which is consistent with the results of previous research [16]. Additionally, the integrated LS-combined AR strategy obtains position solutions more accurately. Specifically, the results of single-frequency PPP obtained for the four frequencies reveal that the accuracy in eastward (E), northward (N), and upward (U) directions is improved by at least 4.1%, 9.6%, and 1.2%, respectively, for B1I; 6.6%, 5.3%, and 0.2%, respectively, for B3I, 12.5%, 14.3%, and 3.8%, respectively, for B1C; and 5.9%, 7.4%, and 18.1%, respectively, for B2a, relative to the results obtained using the traditional method (none). However, because of the high-quality observation of BDS-3, the integrated strategy is slightly more optimized than the two-step method. The ultimate goal of the present study is the development of rapid and high-accuracy positioning services of BDS-3. Therefore, real-time double-frequency BDS-3 PPP is tested with different strategies for multipath errors in the following experiments.The real-time PPP solution is obtained according to the stage of CUM. A flowchart of real-time BDS-3 double-frequency PPP has been published in the literature [2]. Before obtaining the solution of real-time PPP, the satellite ultrarapid orbit and clock offsets are first estimated [2], which is not the focus of the present study. Additionally, the experimental session of real-time PPP is set as DOY 221–228, 2020, where the combinations of B1I+B3I and B1C+B2a are tested for four iGMAS and four MGEX stations. To avoid the restrictions of the real-time stream in the experiments, the RINEX files of observations are input into real-time PPP solutions. Similarly, four schemes of comparative experiments on BDS-3 multipath delays are designed, namely, no correction (none), correction with the AR model (AR), correction with the two-step LS+AR model (LS+AR), and the integrated estimation of the LS+AR model (one step). However, in comparison with BDS-3 single-frequency PPP, the AR strategy is adopted for the estimated coefficients of the previous day and the multipath delays are mitigated by the predicted values. In contrast, the one-step and two-step strategies are based on the prior information of previous model coefficients. Furthermore, the updated model coefficients are estimated with short-term observations (of an hour or so).The position residuals (BJF1) for three continuous days (DOY 226–228, 2020) are plotted in Figure9 (B1I+B3I) and Figure 10 (B1C+B2a) to show the results of BDS-3 real-time PPP. Moreover, the average daily RMS residuals for eight days, four stations, and three strategies (none, two-step, and one-step) and the corresponding improvements are listed in Table 6. Figures 9 and 10 show that there are no appreciable differences in the final precision among the four strategies, although the one-step correction performs slightly better that the two-step strategy. Furthermore, there is a bias in the direction of U, which might be due to the uncalibrated PCV and phase center offset of the satellites.Figure 9
Position residuals of BDS-3 double-frequency real-time PPP forB1I+B3I (BJF1).Figure 10
Position residuals of BDS-3 double-frequency real-time PPP forB1C+B2a (BJF1).Table 6
Average daily RMS values (m) of position residuals for four stations and three multipath error modeling strategies and corresponding improvements.
StationsDirectionsB1I+B3IB1C+B2aNoneTwo-stepImp. (%)One-stepImp. (%)NoneTwo-stepImp. (%)One-stepImp. (%)BJF1E−0.054−0.04025.9−0.03927.8−0.056−0.04814.3−0.04126.8N−0.040−0.02342.5−0.01660.0−0.043−0.03225.6−0.02346.5U−0.174−0.1579.8−0.09445.9−0.156−0.12718.6−0.09241.0KNDYE−0.055−0.0550.0−0.04321.8−0.047−0.0438.5−0.03819.1N−0.040−0.01660.0−0.01367.5−0.040−0.02537.5−0.02147.5U−0.204−0.09354.4−0.04378.9−0.161−0.1497.5−0.10336.0WUH1E−0.043−0.0430.0−0.03713.9−0.039−0.02633.3−0.02633.3N−0.035−0.01654.3−0.01362.9−0.039−0.02243.6−0.01756.4U−0.183−0.10045.4−0.04674.9−0.154−0.13114.9−0.09240.3XIA1E−0.052−0.01473.1−0.01375.0−0.051−0.0501.9−0.03433.3N−0.054−0.01572.2−0.01375.9−0.033−0.02815.2−0.01651.5U−0.108−0.05846.3−0.02874.1−0.142−0.1345.6−0.10823.9There are no obvious differences in the simulated final position residuals between single-frequency PPP and real-time PPP. Table6 shows that the improvements in E, N, and U directions are at least 13.9%, 60.0%, and 45.9%, respectively, and 19.1%, 46.5%, and 23.9%, respectively, for B1I+B3I and B1C+B2a, respectively. It is suggested that the one-step strategy with the prior information for modeling the BDS-3 multipath has better real-time application than the traditional method (none). It should be noted that the convergence time is also an important indicator for the PPP performances. In experiments, the convergence time and position residuals are simultaneously considered. However, it is suggested that there is an average improvement within 2% for E and N directions based on the proposed method, which cannot be as a significant improvement in experiments. However, it is difficult to count the convergence time of the U direction as a systematic error found. Therefore, the convergence time of the PPP solution is ignored in experiments.More BDS-3 observations collected by GNSS receivers are used to test the proposed strategy of mitigating the BDS-3 multipath delay. Furthermore, the results of the SPP experiment are verified using a set of real observations, the trajectory of which is shown in Figure11. In consideration of the decoding ability of GNSS receivers, the code observations of B1C are analyzed with an interval of 0.2 s, where B1I and B3I are used to model the code multipath delays. In experiments, three test scenarios are conducted with the same observations, strategies, and parameter estimations. However, the BDS-3 code multipath delays are modeled and corrected with the proposed method (one-step), traditional method (two-step), and none considerations (none) in the preprocessing of BDS-3 observations. In a word, the changes of position residuals are caused by the different multipath mitigation method. The single-epoch solution of SPP in E and N directions is shown for three strategies (none, one-step, and two-step) in Figure 12. In analyzing the accuracy of SPP, the station positions of GPS double-frequency kinematic PPP are taken as a reference. SPP does not provide the desired result owing to the complex kinematic environments around the stations. However, the precision and stability can be improved easily by adopting the one-step and two-step strategies. Moreover, there are no appreciable differences in results between the one-step and two-step strategies as no prior information is used in modeling BDS-3 multipath errors. Therefore, a more precise and accurate strategy should be developed to outperform the proposed one-step correction for multipath delays in real-time applications. Figure 6 shows that the multiple paths resulting from the kinematic environment cannot be precisely corrected for, whereas the systematic and trend terms along the elevation angle can be reduced using our method. Because of the complex environments, the improvements of dynamic multipath mitigation for position residuals cannot totally be ascribed to the proposed method as the significant errors in different directions. In the next step, the dynamic multipath mitigation will be paid more attention to improve the applications of BDS-3 real-time services.Figure 11
Trajectory of the SPP experiment.Figure 12
Position residuals obtained in SPP experiments with different strategies.


## 4. Discussion

The present study modeled and mitigated the multipath delays of BDS-3 code observations adopting a new strategy, namely, the one-step mitigation strategy, where three points are optimized in contrast with the algorithm traditionally used. First, to avoid the effects of noise and limited observations in modeling the multipath series, sparse modeling is adopted using a combination of TV and EN terms. Second, the multipath series of BDS-3 observations can be divided into two parts, namely, random and systematic components. Therefore, a one-step estimation strategy based on theLS+AR model was proposed. Third, in consideration of real-time and near real-time high-precision applications, the model coefficients of multipath series are estimated and improved with the aid of a prior information.To verify the proposed strategy, four groups of experiments were conducted with one-month iGMAS and MGEX observations. The model residuals suggest that sparse modeling by adding TV and EN terms improves the accuracy of multipath series. Meanwhile, single-frequency PPP and double-frequency PPP were tested for the proposed one-step mitigation method. The accuracy in E, N, and U directions was improved by at least 4.1%, 9.6%, and 1.2%, respectively, for B1I; 6.6%, 5.3%, and 0.2%, respectively, for B3I; 12.5%, 14.3%, and 3.8%, respectively, for B1C; and 5.9%, 7.4%, and 18.1%, respectively, for B2a relative to the use of the traditional (AR) method in BDS-3 single-frequency PPP. Moreover, there were improvements in E, N, and U directions of at least 13.9%, 60.0%, 45.9%, respectively, forB1I+B3I and 19.1%, 46.5% and 23.9%, respectively, for B1C+B2a. Furthermore, the solution of SPP can be optimized using the proposed method. The use of the proposed method thus improves the accuracy of BDS-3 multipath mitigation and real-time application.As discussed above, the multipath delays of BDS-3 code observations affect the PPP solution, especially in single-frequency applications. In general, continuously operating reference stations and base stations can ignore the effects of multipath delays by adopting a suitable cutoff elevation angle. However, three problems relating to the raw observations remain: (1) iGMAS and MGEX code observations suggest that there are multipath delays in different observations, with which values are lower than that of actual situations [20], (2) the multipath delays have systematic and random components, which can be explained by code bias [2] and the complex environments around stations [20], and (3) the most destructive near-field multipath has a frequency exceeding 1 Hz, as verified by real-time SPP solutions. In fact, widely used 30 s observations cannot avoid the effects of the environment and the characteristics of signals. In improving the performance of BDS-3, we need to analyze and mitigate the effects of multiple paths of code observations, which are the result of complex factors, especially for the newly deployed BDS-3 satellites. Furthermore, the method can be used in the application of low-cost receivers [31].Researches about mitigating GPS multipath delays from day-to-day range measurements [32] involved the pseudo-observation with an empirical sample variance, in which systems of errors are divided into two subsystems. Compared with this research, improved algorithms and more analyses are included in our study, in which differences can be summarized as three points: (1) the new-generation BDS-3 satellites and its multifrequency observations are fully discussed. In fact, there are a lot of methods, strategies and algorithms about the mitigation of GPS multipath delays, while the BDS-3 should be carefully studied as its new services, especially for new observations, (2) the purpose of this research mainly focuses on the real-time or near real-time applications, in which the sample mean and variance cannot be used in estimation of multipath parameters. In addition, the BDS-3 code multipath delays are divided into two parts, namely, systematic and random components, to model and reduce the next-arc multipath delays, respectively, and (3) the improved method is more easier to conduct and insert into the preprocessing of PPP software. Moreover, the prior information of unknown coefficients can optimize the estimation models. In general, the new method is further improved and simplified than that of traditional algorithms.In summary, the main goal of the present study was the mitigation of the multipath delay of BDS-3 code observations with the aid of sparse modeling and integrated estimation. We first analyzed the characteristics of BDS-3 multipath delays, which were then modeled using regularization, EN, and TV terms and an integrated LS-combined AR algorithm based on static station observations. Multipath delay coefficients cannot be used for accurate corrections in the real-time and near near-time applications of BDS-3. We therefore proposed an improved method that takes the estimated coefficients as prior information and feeds the coefficients into the real-time correction of BDS-3 raw observations. Because of the function between the multipath delay and elevation, the next-arc raw observations can be further improved with the real-time estimated coefficients through the accumulation of short-arc observations. Similarly, in dynamic observations, we can accumulate static observations over a short period (of an hour or so) to calculate multipath coefficients, while the prior values are set at zero.


## 5. Conclusions

The BDS-3 system adopts new frequencies, new observations, and new technologies in providing global services, and these need to be further discussed for high-precision applications. The performance of the BDS-3 system is important in the research field of GNSS navigation and positioning while innovations and developments of the BDS-3 system are important in further expanding high-quality PNT services. In this paper, the mitigation model of the BDS-3 multipath error was optimized for rapid high-accuracy applications. The PPP solution cannot be improved appreciably because the processing of the multiple frequencies and errors of observations still refers to BDS-2 [2]. In future work, we will consider ultrarapid BDS-3 satellite products and bias parameters in optimizing rapid BDS-3 services.



---

*Source: 9947704-2021-09-23.xml*