Characteristics of BeiDou Navigation Satellite System Multipath and Its Mitigation Method Based on Kalman Filter and Rauch-Tung-Striebel Smoother

Global Navigation Satellite System (GNSS) carrier phase measurement for short baseline meets the requirements of deformation monitoring of large structures. However, the carrier phase multipath effect is the main error source with double difference (DD) processing. There are lots of methods to deal with the multipath errors of Global Position System (GPS) carrier phase data. The BeiDou navigation satellite System (BDS) multipath mitigation is still a research hotspot because the unique constellation design of BDS makes it different to mitigate multipath effects compared to GPS. Multipath error periodically repeats for its strong correlation to geometry of satellites, reflective surface and antenna which is also repetitive. We analyzed the characteristics of orbital periods of BDS satellites which are consistent with multipath repeat periods of corresponding satellites. The results show that the orbital periods and multipath periods for BDS geostationary earth orbit (GEO) and inclined geosynchronous orbit (IGSO) satellites are about one day but the periods of MEO satellites are about seven days. The Kalman filter (KF) and Rauch-Tung-Striebel Smoother (RTSS) was introduced to extract the multipath models from single difference (SD) residuals with traditional sidereal filter (SF). Wavelet filter and Empirical mode decomposition (EMD) were also used to mitigate multipath effects. The experimental results show that the three filters methods all have obvious effect on improvement of baseline accuracy and the performance of KT-RTSS method is slightly better than that of wavelet filter and EMD filter. The baseline vector accuracy on east, north and up (E, N, U) components with KF-RTSS method were improved by 62.8%, 63.6%, 62.5% on day of year 280 and 57.3%, 53.4%, 55.9% on day of year 281, respectively.


Introduction
Recently, the application of Global Navigation Satellite System (GNSS) technology in deformation monitoring has become an important way to monitor the structural health of buildings for its advantages in automation, all-weather, real time and large scale, etc. [1]. Affected by orbit error, ionospheric delay, tropospheric delay, phase center offset and multipath error, the accuracy of GNSS carrier phase measurement generally only reaches centimeter level. To achieve the millimeter or even submillimeter accuracy which is necessary for deformation monitoring, strict elimination of errors is essential. In the case of short baseline (<3 km), most errors excluding multipath can be eliminated to a large extent by the carrier phase double difference (DD) technique because of the high space correlation.
During the last two decades, many GPS multipath mitigation methods have been developed which can be classified to three classes [2]. The first one is selecting an open area with no reflection or shelter for GNSS measurement while it is limited by objective conditions in deformation monitoring. The second one is improving the hardware devices of GNSS receivers to exclude multipath errors when the GNSS signal is used. The last and most common method is data post-processing approaches, including weight adjustment based on data signal to noise ratio (SNR) techniques [3]; ray-tracing approach [4]; and sidereal filtering (SF) based on the coordinate domain [5,6]/observation domain [7,8] method, which is the most popular.
At present, the BeiDou Navigation Satellite System (BDS) consisting of geostationary earth orbit (GEO) satellites, inclined geosynchronous orbit (IGSO) satellites and medium Earth orbit (MEO) satellites can provide high-accuracy location services, which are comparable to Global Position System (GPS) in the Asian-Pacific region [9,10]. Many literatures have researched the performance of the combination of BDS/GPS [11,12]. By 2020, the service range of BDS will cover all the world and its data processing has become a hot spot and focus of GNSS positioning. It is necessary to conduct accuracy assessments and to develop error cancellation techniques about BDS positioning such as multipath mitigation.
Since the relative geometry of a GPS satellite with receiver and reflective surface repeats itself in each sidereal day, multipath errors are highly correlated between consecutive sidereal days if the geometry changes are the same [13,14]. Hence, the SF based on coordinate domain using daily repetition of multipath can extract the multipath model from the coordinate sequence of the first day to correct that of the next day at a fixed station. However, this method uses the average satellite orbital period and ignores the differences between satellites, so it cannot be effectively applied in BDS multipath mitigation with its three different types of orbital satellites. Axelrad et al. found that the orbital repeat periods of GPS satellites are not strictly a sidereal day (about 23 h and 56 min) and there is a time shift which is not a constant for different satellites [15]. Therefore, the SF based on observation domain, which calculates the orbital period of each satellite, extracts and corrects the multipath model on corresponding satellite observations between the observation sequences of two multipath periods, was proposed and developed. Zhong et al. developed a SF based on single difference (SD) residuals for mitigating GPS multipath effects on short baselines [7]. Ye et al. studied the multipath repeat cycle of BDS satellites, and proposed a SF carrier phase multipath elimination approach of BDS system [14]. Dong et al. proposed a Multipath Hemispherical Map (MHM) to achieve real-time resolution and correction of multipath errors [16]. Dai et al. compared the MHM algorithm with the sidereal filtering algorithm, and proposed an improved multipath error parameterization model for BDS GEO multipath elimination [17].
In order to obtain the multipath error model, it is essential to select an effective filtering method which controls the effect of multipath elimination. Lau proposed a robust three-level wavelet packet for carrier phase multipath [18]. Zhong et al. proposed an adaptive wavelet transform based on a cross-validation method to mitigate GPS multipath effects [19]. Dai et al. applied empirical mode decomposition (EMD) in denoising coordinate sequence of short GPS baseline [20]. Kaloop et al. used four prediction models which are applied and used with neural network solutions: back-propagation, Cascade-forward back-propagation, adaptive filter and extended Kalman filter, to extract GPS time-series model and denoise the data for monitoring of structures [21]. FIR filtering has also been used in extracting multipath models in literature [22]. To obtain a more accurate multipath model, a multipath model extraction method based on Kalman filter and Rauch-Tung-Striebel Smoother (RTSS) was proposed and the multipath repetition characteristics of BDS satellites were explained by qualitative analysis and quantitative calculation in the paper. To verify the performance of the new method, we compared it with the traditional wavelet and EMD methods used in [19,20].
The paper is arranged as follows: in Section 2, the basic principle of multipath error is explained and the steps of reconstructing single difference (SD) residuals is shown with a flowchart. Qualitative analysis and quantitative calculation of orbital periods of BDS satellites are carried out in Section 3. Section 4 introduces three filter methods for denoising SD residuals for extracting multipath models and a simulation experiment was carried out. In Section 5, we verified the consistency between orbital periods and multipath periods of BDS satellites through a set of measured data, applied three filtering methods in extracting multipath models and found they all achieve good results in multipath mitigation. Finally, conclusions were drawn in Section 6.

Basic Principle
The multipath problem is a widespread problem in signal propagation, such as the network signal multipath and the sensor data delivery multipath including the underwater sonar multipath and the satellite signal multipath. Data delivery path in Information-Centric Sensor Networks was studied in [23,24]. The AUV Positioning Method Based on Tightly Coupled SINS/LBL was proposed to deal with the Underwater Acoustic Propagation Multipath in [25]. In [26], the Kalman filter was introduced to realize multi-path network synchronization. Satellite signal multipath mitigation is a key step in precise GNSS positioning and is our focus work. The path deviation of satellite signals caused by multipath effects can be expressed as [7]: where D represents the horizontal distance between the receiver antenna and the reflector around the station, θ and λ are respectively the incidence angle and wavelength of the reflected signal, α is the reflection coefficient of reflector and ∆ϕ m represents the phase delay caused by multipath effects. It can be found that the multipath errors are mostly affected by horizontal distance D, incident angle θ and reflection coefficient α. However, the reflective surface, antenna and their properties remain unchanged in case of a static baseline, which means that the horizontal distance and reflection coefficient are constant. Hence it is only the incidence angle of reflected signals that could change multipath errors. Meanwhile, the incidence angle of reflected signals is determined by position of satellites in the sky which would periodically repeat in specific orbit, so multipath errors are also periodic. It indicates that we can mitigate multipath effects for current data by subtracting the multipath error models extracted from the last period data.

Integer Ambiguity Estimation and Validation with Multipath Errors
The integer ambiguity estimation and validation are the two basic steps of high precision carrier phase positioning. The least-squares ambiguity decorrelation adjustment (LAMBDA) has always been the most popular algorithm for ambiguity estimation since Teuniseen proposed it in 1990s [27]. The core idea is to reduce the search space by decorrelating the covariance matrix of float ambiguity to improve the search efficiency. Many modified algorithms were developed based on LAMBDA, such as the mixed upper/lower triangular integer Cholesky decomposition algorithm [28], inverse integer Cholesky decomposition algorithm [29] PAR LAMBDA algorithm [30]. Another representative algorithm called MLAMBDA algorithm was proposed by Chang [31] and is used for ambiguity resolution in this paper.
Another key step is ambiguity validation which can ensure the reliability of fixed ambiguity. The common methods include three types: ratio test based on statistics [32], index method of success rate/failure rate [33] and the combination of the two previous methods [34]. In this paper, the R-ratio test was applied to validate the ambiguity. The DD equations without/with multipath errors can be expressed as: Sensors 2018, 18, 198 4 of 20 where y/y' represent the GNSS observations without/with multipath errors, a and b represent the ambiguity vector and baseline vector with coefficient matrixes A and B, m and e are the multipath errors and DD residuals, respectively. Classifying multipath errors into DD residuals, we get: where e' = e − m, then −m is the DD multipath errors remaining in DD residuals e'. We can compute the float solutions, integer ambiguity and fixed solution using the next criteria: where ||M|| 2 N = M T N M and Q −1 y and Q −1 a are the variance-covariance matrixes of DD observations and DD float ambiguity, respectively. The main contribution of the LAMBDA algorithm is introducing a Z-transformation [27] by z = Z T a and Qẑ = Z T QâZ, then criterion (6) is changed to: Here ambiguity correlation is reduced and search efficiency is improved. After integer ambiguity estimation the ratio test can be carried out with the following criterion [32]: where and is the quadratic form of ambiguity residuals, Ω min and Ω sec represent the minimum and secondary values of Ω i , respectively. The threshold is set to 3 and if ratio ≥ 3, the integer ambiguityȃ corresponding to Ω min is considered to be fixed correctly [34].
Actually, the maximum effect on carrier phase by multipath errors is only up to 1/4 wavelengths when arctan( α sin( 4πD λ sin θ) 1+α cos( 4πD λ sin θ) ) approaches π 2 according to Equation (1). Two cases should be considered in ambiguity estimation: first, multipath errors are so small that they cannot change the estimation of integer ambiguity-but could affect the result of the baseline vector by changing phase observations (the non-integer part of the phase); second, multipath errors bring out a carrier phase measurement error of a few centimeters and both the ambiguity estimation and baseline vector are changed. Therefore, both two cases will be improved if multipath errors were removed from residuals.

Reconstruction of SD Residuals and Correction of Multipath Error
Since most model errors excluding multipath error can be eliminated by DD processing, multipath error becomes the main component of the remaining errors in DD observations. There are two disadvantages when DD residuals are applied to obtain multipath models. Firstly, the reference satellite in the DD process is the satellite with the highest elevation angle, which means that once the reference satellite changes, the DD residuals sequence will become complex and cannot correspond to the specific satellites. Secondly, it is difficult to determine the orbital period of two satellites as their orbital periods are different from each other. Reference [35] gives another way that we can reconstruct SD residuals from DD residuals by the transition matrix which contains the "zero mean" assumption ∑ w i s i AB = 0, then the reconstructed SD residuals can be expressed as: where s i AB and d ij AB represent reconstructed SD residuals and DD residuals, respectively; w i is the weighting factor. Because we use the weighting model of elevation angles, w i = sin 2 (θ) here. Then the process of calculating SD residuals and correcting multipath errors is as shown in Figure 1: (1) The mean solution of the baseline vector of first period data was calculated as the true baseline vector b; (2) The baseline solution is substituted back to the DD equations of first period data to compute ambiguity inversely. Then the float ambiguity solution is fixed to integer ambiguity a; (3) Next, both b and a were substituted into the DD equations to obtain DD residuals consisting of multipath errors and Gaussian white noise; (4) DD residuals were reconstructed to SD residuals by Equation (10); (5) The multipath model of each BDS satellite was extracted from SD residuals through different filter methods with random noise being filtered; (6) The multipath model sequences were adapted according to the time shift computed from first period data; (7) The adaptive multipath model sequences were substituted into the SD equations of the second period data to mitigate multipath effects according to high correlation between multipath errors of two period data; (8) Finally, the baseline coordinate sequence of the second period data was evaluated from corrected DD equations.
where and represent reconstructed SD residuals and DD residuals, respectively; wi is the weighting factor. Because we use the weighting model of elevation angles, wi = sin 2 (θ) here. Then the process of calculating SD residuals and correcting multipath errors is as shown in Figure 1: (1) The mean solution of the baseline vector of first period data was calculated as the true baseline vector b; (2) The baseline solution is substituted back to the DD equations of first period data to compute ambiguity inversely. Then the float ambiguity solution is fixed to integer ambiguity a; (3) Next, both b and a were substituted into the DD equations to obtain DD residuals consisting of multipath errors and Gaussian white noise; (4) DD residuals were reconstructed to SD residuals by Equation (10); (5) The multipath model of each BDS satellite was extracted from SD residuals through different filter methods with random noise being filtered; (6) The multipath model sequences were adapted according to the time shift computed from first period data; (7) The adaptive multipath model sequences were substituted into the SD equations of the second period data to mitigate multipath effects according to high correlation between multipath errors of two period data; (8) Finally, the baseline coordinate sequence of the second period data was evaluated from corrected DD equations.

Qualitative Analysis
All the GPS satellites are the same orbital satellite-MEO satellite and their orbital periods are almost one sidereal day although there may be a deviation of several minutes (notice that the orbital period here is not calculated in the celestial coordinate system, but in the site-centric coordinate system) However, for the three types of satellites of BDS system, the orbital repeatability is more complex. Figure 2 shows the satellite tracks of BDS GEO satellite (C01), IGSO satellite (C06) and MEO

Qualitative Analysis
All the GPS satellites are the same orbital satellite-MEO satellite and their orbital periods are almost one sidereal day although there may be a deviation of several minutes (notice that the orbital period here is not calculated in the celestial coordinate system, but in the site-centric coordinate system) However, for the three types of satellites of BDS system, the orbital repeatability is more complex. Figure 2 shows the satellite tracks of BDS GEO satellite (C01), IGSO satellite (C06) and MEO satellite (C11) on day of year (DOY) 284, 285 and 291 in 2016. We can find that the ground tracks of the C06 is like an "8" shape and its tracks on the three days almost coincide. This indicates that the satellite orbital periods of IGSO satellites are about 1 day. Meanwhile, the tracks of C11 just coincide on DOY 284 and 291, but that on DOY 285 is independent of the other two days. Consequently, the orbital period of BDS MEO satellites is about 7 days. In Figure 2, we can also find that C01 almost remained stationary during the three days, but actually it still moved in a very small area because it was vulnerable to tidal changes, atmospheric pressure and nutation and other external intrusion. As shown in Figure 3, BDS GEO satellites will move in the shape of a small "8" within 4 degrees in latitude and 1 degree in longitude. It is clear that the orbital periods of GEO satellites are about 1 sidereal day. We can find that the ground tracks of the C06 is like an "8" shape and its tracks on the three days almost coincide. This indicates that the satellite orbital periods of IGSO satellites are about 1 day. Meanwhile, the tracks of C11 just coincide on DOY 284 and 291, but that on DOY 285 is independent of the other two days. Consequently, the orbital period of BDS MEO satellites is about 7 days. In Figure 2, we can also find that C01 almost remained stationary during the three days, but actually it still moved in a very small area because it was vulnerable to tidal changes, atmospheric pressure and nutation and other external intrusion. As shown in Figure 3, BDS GEO satellites will move in the shape of a small "8" within 4 degrees in latitude and 1 degree in longitude. It is clear that the orbital periods of GEO satellites are about 1 sidereal day.   We can find that the ground tracks of the C06 is like an "8" shape and its tracks on the three days almost coincide. This indicates that the satellite orbital periods of IGSO satellites are about 1 day. Meanwhile, the tracks of C11 just coincide on DOY 284 and 291, but that on DOY 285 is independent of the other two days. Consequently, the orbital period of BDS MEO satellites is about 7 days. In Figure 2, we can also find that C01 almost remained stationary during the three days, but actually it still moved in a very small area because it was vulnerable to tidal changes, atmospheric pressure and nutation and other external intrusion. As shown in Figure 3, BDS GEO satellites will move in the shape of a small "8" within 4 degrees in latitude and 1 degree in longitude. It is clear that the orbital periods of GEO satellites are about 1 sidereal day.

Quantitative Calculation
After qualitative analysis of the BDS satellites orbital periods, quantitative calculation for each satellite is also needed. On one side, the orbital periods of satellites are not strictly integer days and there is a time shift. On the other hand, the orbital periods of BDS satellites are different from each other. According to the orbit radius and the angular velocity correction given by the broadcast ephemeris, the angular velocity of the satellite can be obtained at each epoch: where GM = 1.996498 × 10 7 is the square root of the product of the gravitational constant of the earth and the mass of the earth, ∆n is the correction parameter of satellite angular velocity, a is the square root of the long half axis of the satellite orbit ellipse. For BDS GEO and IGSO satellites, the time shifts can be computed by: We also obtain the time shifts of MEO satellites by: Based on the above formulas and BDS ephemeris in January and February 2017, we calculated the time shifts of all satellites at each ephemeris epoch as shown in Figure 4 and Table 1. From Figure 4, it is easy to find that the time shifts of the three different types of satellites are different. The corresponding time shifts of different satellites with the same type varies at different time. Even the time shift of one satellite varies at different time. As Figure 4a shows, the time shift curves of GEO satellites fluctuate greatly even in one day and the range is about 30 s. The result also indicates that GEO satellite has a high-frequency jitter during operation due to interference from external forces. Moreover, it is easy to find that the time shifts of GEO1-GEO5 satellites come into a similar variation trend in the two months so their mean value are nearly the same values with 235 s or 236 s. In Figure 4b we can find that, the time shifts of the IGSO satellites is relatively stable in one day and there is no great volatility like GEO satellites. However, the curves of time shifts of different IGSO satellites keep independent of each other, consequently, a few seconds of deviation occurs among the time shifts of IGSO satellites. The minimum value is 233 s corresponding to IGSO9 and the maximum value is 245 s corresponding to IGSO10. Only IGSO8 happens to a great jump about 50 s on 22 January. The main reason for this phenomenon may be orbital maneuver. For the MEO satellites in Figure 4c, the change tendency of each curve fluctuates is around 1700 s although there still exist differences. However, it can be found from Table 1 that the range of time shift of MEO satellites are 12-16 s, which is smaller than that of GEO and IGSO satellites.

Quantitative Calculation
After qualitative analysis of the BDS satellites orbital periods, quantitative calculation for each satellite is also needed. On one side, the orbital periods of satellites are not strictly integer days and there is a time shift. On the other hand, the orbital periods of BDS satellites are different from each other. According to the orbit radius and the angular velocity correction given by the broadcast ephemeris, the angular velocity of the satellite can be obtained at each epoch: where GM = 1.996498  10 7 is the square root of the product of the gravitational constant of the earth and the mass of the earth, n  is the correction parameter of satellite angular velocity, a is the square root of the long half axis of the satellite orbit ellipse. For BDS GEO and IGSO satellites, the time shifts can be computed by: We also obtain the time shifts of MEO satellites by: Based on the above formulas and BDS ephemeris in January and February 2017, we calculated the time shifts of all satellites at each ephemeris epoch as shown in Figure 4 and Table 1. From Figure  4, it is easy to find that the time shifts of the three different types of satellites are different. The corresponding time shifts of different satellites with the same type varies at different time. Even the time shift of one satellite varies at different time. As Figure 4a shows, the time shift curves of GEO satellites fluctuate greatly even in one day and the range is about 30 s. The result also indicates that GEO satellite has a high-frequency jitter during operation due to interference from external forces. Moreover, it is easy to find that the time shifts of GEO1-GEO5 satellites come into a similar variation trend in the two months so their mean value are nearly the same values with 235 s or 236 s. In Figure  4b we can find that, the time shifts of the IGSO satellites is relatively stable in one day and there is no great volatility like GEO satellites. However, the curves of time shifts of different IGSO satellites keep independent of each other, consequently, a few seconds of deviation occurs among the time shifts of IGSO satellites. The minimum value is 233 s corresponding to IGSO9 and the maximum value is 245s corresponding to IGSO10. Only IGSO8 happens to a great jump about 50 s on 22 January. The main reason for this phenomenon may be orbital maneuver. For the MEO satellites in Figure 4c, the change tendency of each curve fluctuates is around 1700 s although there still exist differences. However, it can be found from Table 1 that the range of time shift of MEO satellites are 12-16 s, which is smaller than that of GEO and IGSO satellites.

Filter Methods for Extracting Multipath Error Models
As mentioned in Section 2.3, the DD and SD residuals consist of multipath errors and white noise. The indispensable step of multipath mitigation is extracting multipath models from SD residuals. Denoising SD residuals to approach multipath errors is the main part of this process where reasonable filter should be selected. Many researchers have studied on filter methods for multipath mitigation such as wavelet filter, EMD filter and FIR filter. In this paper, Kalman filter and RTS smoother (KF-RTSS) was firstly used in BDS multipath mitigation. The basic principle of the wavelet filter and EMD filter were explained briefly and the process of KF-RTSS was described in detail as follows.

Filter Methods for Extracting Multipath Error Models
As mentioned in Section 2.3, the DD and SD residuals consist of multipath errors and white noise. The indispensable step of multipath mitigation is extracting multipath models from SD residuals. Denoising SD residuals to approach multipath errors is the main part of this process where reasonable filter should be selected. Many researchers have studied on filter methods for multipath mitigation such as wavelet filter, EMD filter and FIR filter. In this paper, Kalman filter and RTS smoother (KF-RTSS) was firstly used in BDS multipath mitigation. The basic principle of the wavelet filter and EMD filter were explained briefly and the process of KF-RTSS was described in detail as follows.

Wavelet Filter
A wavelet filter decomposes a signal into different frequency layers by Fourier transform [36]. In most cases, the low frequency layers are the useful part corresponding to the multipath errors of the SD residuals, and the high frequency layers contains white noise that is going to be filtered. Some details of the high frequency layers are still preserved by different threshold methods. Then the preserved details and low frequency part are reconstructed to approach original data. The process of wavelet filter consists of selecting wavelet basis function, determining decomposition level and defining threshold, respectively. The wavelet basis function we used is sym6, the decomposition progression is 4 and the threshold of filter is default provided by MATLAB.

EMD Filter
Empirical mode decomposition (EMD) is proposed by Huang in 1998 [20,37]. As a new signal decomposition method, it is a part of the Hilbert Huang transform (HHT) method which can analyze the spectrum of non-linear and non-stationary signals adaptively. EMD is mainly applied to obtain the instantaneous frequency characteristics of signal sequences. Hence it is adaptive to decompose data into a set of Intrinsic Mode Function (IMF) components with different frequencies, which represent the basis of the data. Similar to the wavelet filter, the high-frequency IMF component is considered as the noise part. Although the wavelet filter has multi-resolution and high time-frequency domain resolution, its application is limited by many problems, such as wavelet basic selection and decomposition progression determination. Reference [38] shows that the product of the energy density and the mean period of the IMF component of the Gaussian white noise signal is a constant. Hence we can determine the noise progression by Equation (14): where ET k is the product of the energy density and the mean period of the IMF k, R k is the ratio factor that control noise progression. Generally, when R k ≥ 2, the sum of the first IMF to IMF k − 1 is the noise part of the original signal and the sum of other IMF components is the filtered signal. Therefore, EMD offers a different approach to data decomposition, and we can apply it to the study of extracting BDS multipath models.

Kalman Filter and RTS Smoother
Kalman filter is the traditional algorithm for the Kalman filter (KF) and Rauch-Tung-Strieber smoother (RTSS) can provide the minimum mean square estimates (MMSE) of states for state-space models with additive Gaussian system and observation noises given a series of the past, current, and future observations [39]. Hence, it has been widely used in data denoising processing. Wang applied KF-RTSS method in the position estimation of a GPS software receiver [40]. Sarkka [42]. In this paper, the KF-RTSS method is introduced to extract the BDS multipath models from SD residuals and correct multipath errors. The simplest state space model is a linear model, which can be expressed as follows: where x k is the system state of k step, z k is the measurement of k step, F k−1 is the transition matrix transforming the system state from the k − 1 step to the K step and H k is the transition matrix giving the relation between the observation vector z k and the state vector x k . w k−1 and v k−1 are the process noise and measurement noise of state equation and measurement equation, respectively. Meanwhile, these are theoretically white noise that meets ε(w k−1 ) = 0, ε(v k ) = 0 and ε w k v T k = 0, ε w 0 w T k = 0, ε x 0 v T k = 0, which means measurement noise is irrelevant to system status, process noise is irrelevant to observation, and measurement noise is also irrelevant to process noise. Let ε w k−1 w T k−1 = Q k−1 and ε v k v T k = R k , then the iterative process of KF and RTSS is as follows: (1) Prediction and update of state and measurement (KF) One-step prediction:x Variance matrix of one-step prediction: The state estimation is corrected by the measurement error term: Update the variance matrix of state: where K k represents the filter gain matrix of step k, which can be expressed as: (2) After step k, the state estimator of the K step dynamic system is smoothed by the following formulas (RTSS):x where,x k+1 andP k+1 are the smoother estimates for the state mean and state covariance on time step k. x S k and P S k are the filter estimates for the state mean and state covariance on time step k.x k+1 and P k+1 are the predicted state mean and state covariance on time step k + 1, which are the same as in the Kalman filter. C k is the smoother gain on time step k, which tells how much the smoothed estimates should be corrected on that particular time step.
The process of KF and RTSS is shown as Figure 5. The input SD residuals are the measurement z k . The filtered and smoothed output of multipath models are the state x S k . The system state is filtered from step 1 to k and then each step is smoothed conversely according to filtered value and smoothed value of next step. The initial value of the smoothing process x S k is equal to the last state estimation of the filtering process x k .
S k Figure 5. Process of KF-RTSS method.

Simulation Experiment
Because the multipath error sequence contains signals with different frequencies, we can simulate multipath signal by the combination of three sine wave signals A Gauss white noise with a standard deviation of 1 is added, then the data model turns to: where, Mt is the mixed signal which can be considered as BDS SD residual sequence, et is the Gauss white noise which should be filtered. The sampling interval is 1 s, and the sampling number is 5000. In order to verify the denoising effect of KF-RTSS method, wavelet filter and EMD filter are also used as the comparison. As shown on the left of Figure 6, from top to bottom are the original signal, Gauss white noise and noisy signal. On the right, extracted signals by wavelet filter, EMD filter and KF-RTSS method are shown respectively. Table 2 gives the correlation coefficient between filtered signals and the original signal. The root mean square errors of both signals are assessed. It is clearly that the correlation coefficient of KF-RTSS signal and original signal is the highest (0.9927) while RMSE of EMD signal is closest to original signal. In the three filters, Wavelet signal show the worst performance. Generally, correlation coefficients between filtered signals and original signal are all more than 0.97, which indicates that three filters all have good performance on denoising the contaminated signals.

Simulation Experiment
Because the multipath error sequence contains signals with different frequencies, we can simulate multipath signal by the combination of three sine wave signals S t = sin(2πt/200) + sin(2πt/400) + sin(2πt/600) (28) A Gauss white noise with a standard deviation of 1 is added, then the data model turns to: where, M t is the mixed signal which can be considered as BDS SD residual sequence, et is the Gauss white noise which should be filtered. The sampling interval is 1 s, and the sampling number is 5000. In order to verify the denoising effect of KF-RTSS method, wavelet filter and EMD filter are also used as the comparison. As shown on the left of Figure 6, from top to bottom are the original signal, Gauss white noise and noisy signal. On the right, extracted signals by wavelet filter, EMD filter and KF-RTSS method are shown respectively. Table 2 gives the correlation coefficient between filtered signals and the original signal. The root mean square errors of both signals are assessed. It is clearly that the correlation coefficient of KF-RTSS signal and original signal is the highest (0.9927) while RMSE of EMD signal is closest to original signal. In the three filters, Wavelet signal show the worst performance. Generally, correlation coefficients between filtered signals and original signal are all more than 0.97, which indicates that three filters all have good performance on denoising the contaminated signals.

Simulation Experiment
Because the multipath error sequence contains signals with different frequencies, we can simulate multipath signal by the combination of three sine wave signals sin(2 / 200) sin(2 / 400) sin(2 / 600) A Gauss white noise with a standard deviation of 1 is added, then the data model turns to: where, Mt is the mixed signal which can be considered as BDS SD residual sequence, et is the Gauss white noise which should be filtered. The sampling interval is 1 s, and the sampling number is 5000. In order to verify the denoising effect of KF-RTSS method, wavelet filter and EMD filter are also used as the comparison. As shown on the left of Figure 6, from top to bottom are the original signal, Gauss white noise and noisy signal. On the right, extracted signals by wavelet filter, EMD filter and KF-RTSS method are shown respectively. Table 2 gives the correlation coefficient between filtered signals and the original signal. The root mean square errors of both signals are assessed. It is clearly that the correlation coefficient of KF-RTSS signal and original signal is the highest (0.9927) while RMSE of EMD signal is closest to original signal. In the three filters, Wavelet signal show the worst performance. Generally, correlation coefficients between filtered signals and original signal are all more than 0.97, which indicates that three filters all have good performance on denoising the contaminated signals.

Performance Analysis of BDS Multipath Mitigation Based on the KF-RTSS Method
In order to verify the consistency between orbital period and multipath periods of BDS satellites and the feasibility of the KF-RTSS proposed, we calculated the SD residuals and extract the multipath model with respect to each satellite, and evaluate the baseline results after three filtering methods. This set of data were collected from 27 September to 10 October 2017 year with the frequency of 1 Hz and 5 • elevation mask angle at the roof of School of environment and Geometrics, China University of Mining and Technology. The two GNSS receivers are Trimble R10 (Trimble Navigation, Sunnyvale, California, United States) units and one of them was set in an unobstructed environment as base station while another one was placed about 5 m away from the southeast direction of a white wall as the rover station. The length of the baseline is about 62.210 m. To ensure that the base station is not affected by multipath effects and the rover station has a clear multipath effect, we turned off the anti-multipath function of the rover station while switching on the function of the base station. The processing of extracting SD residuals and mitigating multipath errors are implemented by the GNSS data processing software developed by ourselves.

SD Residual of Different Types of BDS Satellites
According to the steps in Section 2.3, we calculated the SD residuals of different type of BDS satellites.

Performance Analysis of BDS Multipath Mitigation Based on the KF-RTSS Method
In order to verify the consistency between orbital period and multipath periods of BDS satellites and the feasibility of the KF-RTSS proposed, we calculated the SD residuals and extract the multipath model with respect to each satellite, and evaluate the baseline results after three filtering methods. This set of data were collected from 27 September to 10 October 2017 year with the frequency of 1 Hz and 5° elevation mask angle at the roof of School of environment and Geometrics, China University of Mining and Technology. The two GNSS receivers are Trimble R10 (Trimble Navigation, Sunnyvale, California, United States) units and one of them was set in an unobstructed environment as base station while another one was placed about 5 m away from the southeast direction of a white wall as the rover station. The length of the baseline is about 62.210 m. To ensure that the base station is not affected by multipath effects and the rover station has a clear multipath effect, we turned off the anti-multipath function of the rover station while switching on the function of the base station. The processing of extracting SD residuals and mitigating multipath errors are implemented by the GNSS data processing software developed by ourselves.

SD Residual of Different Types of BDS Satellites
According to the steps in Section 2.3, we calculated the SD residuals of different type of BDS satellites. Figure

Extracting Multipath Error Models from SD Difference Residuals
The SD residuals of the last period data cannot be subtracted from SD observations of current period data because the SD residuals consist of multipath error models and white errors which are interference terms and should be filtered. Multipath models of IGSO6 for DOY 272-274 extracted from SD residuals by the three filters mentioned are shown in Figure 9. Table 5 gives the correlation coefficients of multipath models between DOY 273 and 274 for BDS PRN 1-10 satellites with three filtering methods. Compared with Table 4, we can find that the mean correlation coefficients with wavelet filter, KT-RTSS and EMD filtering are 0.9036, 0.9214 and 0.9158, respectively. These performances are all better than that of SD residuals with the value of 0.853. This indicates that three kinds of filter are all effective, therein the performance of KT-RTSS method is the best.

Extracting Multipath Error Models from SD Difference Residuals
The SD residuals of the last period data cannot be subtracted from SD observations of current period data because the SD residuals consist of multipath error models and white errors which are interference terms and should be filtered. Multipath models of IGSO6 for DOY 272-274 extracted from SD residuals by the three filters mentioned are shown in Figure 9. Table 5 gives the correlation coefficients of multipath models between DOY 273 and 274 for BDS PRN 1-10 satellites with three filtering methods. Compared with Table 4, we can find that the mean correlation coefficients with wavelet filter, KT-RTSS and EMD filtering are 0.9036, 0.9214 and 0.9158, respectively. These performances are all better than that of SD residuals with the value of 0.853. This indicates that three kinds of filter are all effective, therein the performance of KT-RTSS method is the best.

Improvement of Baseline Vector with Different Multipath Mitigation Methods
The last step of multipath mitigation is subtracting the multipath models of last period data from current period data and calculate the baseline vector. Shown in Figure 10 are the BDS baseline vectors before and after multipath mitigation with different methods. Table 6 gives the RMSE of deviations in E, N, U components and percentage of improvement respectively. The percentage of improvement with KF-RTSS method is about 62.8%, 63.6%, 62.5% and 57.3%, 53.4%, 55.9% in E, N, U components on DOY 280 and 281, which is slightly better than that of Wavelet filter and EMD filter.

Conclusions
In order to meet the requirements of deformation monitoring with BDS precise positioning techniques for large structures, it is of significance to mitigate multipath effects in carrier phase observations. In this paper, detailed analysis was provided about the characteristics of orbital periods and multipath repeat periods of BDS satellites. It was found that they are consistent. The period of orbit and multipath errors for BDS GEO and IGSO satellites is about one day and that of BDS MEO satellites is about seven days. Sidereal filters based on SD difference were applied to correct multipath errors. The KF-RTSS method was introduced to extract multipath models from single difference (SD) residuals. After subtracting multipath models, the improvement of baseline accuracy on E, N, U directions are about 62.8%, 63.6%, 62.5% on DOY 280 and 57.3%, 53.4%, 55.9% on DOY 281, respectively. Wavelet filter and EMD filter were also used in multipath mitigation. The experimental results indicate that the three filters all have obvious effect on improvement of baseline accuracy and the performance of KT-RTSS method is slightly better than that of the wavelet filter and EMD filter.