A Multipath Error Reduction Method for BDS Using Tikhonov Regularization with Parameter Optimization

Abstract

Multipath error is an important factor restricting the relative positioning accuracy of the Beidou Navigation Satellite System (BDS). Because of the complexity of the reflection environment, the mathematical modeling of multipath errors is quite difficult. The sidereal filtering algorithm corrects multipath errors by using the feature of period repetition, which can greatly reduce its influence and improve the accuracy of positioning and attitude measurement. In view of the constellation heterogeneity of BDS, it is more complicated to apply sideral filtering. Based on the reconstructed single-difference residual of the carrier phase, the multipath repetition time of the Beidou satellite is estimated using the idea of segmentation. The Tikhonov regularization method and the classical wavelet method are used to extract the multipath of the single-difference residual of the carrier phase, and the “clean” sequence of the single-difference residual is obtained. The experimental results show that it is feasible to extract the multipath error correctly by Tikhonov regularization, and the multipath error is smoother than the original residual measurement. Furthermore, the estimation method of the regularization parameter is further optimized. After using the optimized Tikhonov regularization method with sidereal filtering, the mean RMS improvements of GEO, IGSO, and MEO satellites are 45.9%, 38.2%, and 37.5%, respectively. The positioning accuracy on E, N, and U components is improved by 24.8%, 26.3%, and 42.7%, respectively. The attitude resolution accuracy is improved by 22.9% in the yaw angle and 12.6% in the pitch angle. The proposed method can be an alternative BDS multipath error modeling and mitigation approach.

1. Introduction

Although new positioning techniques and modes are constantly emerging, carrier-relative positioning is still commonly used in high-precision GNSS applications. Multipath error is a crucial error source for high-precision positioning in GNSS. It can be mitigated by selecting appropriate observation station positions, improving receiver hardware equipment, and using data post-processing modeling methods [1,2,3,4,5,6]. The data post-processing modeling method mainly uses information such as the data recorded by the receiver, the reflection environment, the receiver antenna, and the number of visible satellites. Multipath error is extracted from the original observations or results files to weaken its influence on the positioning. It can be summarized into time-domain and spatial-domain modeling methods. The time-domain modeling method takes advantage of the repeatability of the multipath period to weaken the influence of the multipath effect, also known as sideral filtering. Spatial-domain modeling methods include the grid method [7,8], a spherical harmonic method [9,10], the support vector machine method [11,12], etc.

As one of the most common methods of time-domain modeling, sidereal filtering has two key steps in its application. Firstly, it is necessary to establish an accurate multipath extraction model, and widely used techniques include empirical mode decomposition (EMD) [13,14], wavelet filtering [15,16,17], Vondrak filtering [18,19], Kalman filtering [20,21,22], etc. Su Mingkun et al. used adaptive thresholding wavelet denoising and the single difference sideral filtering method with the double reference shift strategy (ATDR) to eliminate the multipath error in the static short baseline. Compared with the constant thresholding method, positioning accuracy is significantly improved by applying the adaptive thresholding method [23]. Zhong Ping et al. used the sideral filtering algorithm of single-difference residuals to weaken GPS multipath error. They verified through experiments that the results based on the observation domain were better than the multipath error correction results based on the coordinate domain [24]. Chang Guobin et al. regard multipath error modeling as a signal reconstruction problem and use Tikhonov first-order difference regularization and second-order difference regularization to extract multipath. Through experiments, it was verified that the first-order and second-order regularization methods have basically the same performance, and both are slightly better than the wavelet denoising method [25]. Secondly, it is very important to determine the multipath repeat time (MRT). The common MRT estimation methods include the orbit repeat time method (ORTM), aspect repeat time adjustment (ARTA), and residual correlation method (RCM) [26,27,28]. The results of the three methods have been compared and analyzed in the existing literature. The effect of RCM directly applied to sideral filtering is slightly better than the other two methods, but the calculation amount of ORTM is minimal [29]. Since satellite orbit is subject to jitter, the repetition period of different GPS satellites varies slightly at different times. Yin Haitao et al. proposed a modified sidereal filtering algorithm based on the segmentation idea, and the positioning accuracy of the modified coordinate components (E, N, and U) was slightly improved [30].

Because of the characteristics of BDS satellite design and constellation composition, the multipath effect characteristics of BDS observations are significantly different from those of other GNSS systems. The effect of GEO and IGSO on BDS multipath error on the final results for positioning and attitude determination remains obvious [31]. In addition, GEO is not strictly stationary but slowly changing, with roughly the same trend for two consecutive days. The IGSO error is similar to GEO, but the magnitude change is larger than GEO. There is no continuity of multipath errors in MEO for two consecutive days [32]. In the case of a fixed observation environment, the multipath error caused by satellite signal reflection is similar to the orbital repetition period of various types of satellites. In other words, the multipath error repetition period of GEO and IGSO is about 1 day, and that of MEO is about 7 days [33].

With this background, this study contributes a Tikhonov regularization method with parameter optimization to extract the multipath error sequence of BDS. It combines the sidereal filtering method to correct the multipath error of each BDS satellite. The correctness and effectiveness of this method are verified by experiments. The rest of this study is organized as follows: Section 2 introduces the theoretical methods and models, including the reconstruction of the single-difference residual of the carrier phase, Tikhonov regularization modeling, optimization of regularization parameters based on the bootstrap method, wavelet filtering, and attitude determination theory. Section 3 shows the details of the experiment and the analysis of the results. Section 4 discusses the proposed method, and Section 5 is our final conclusion to the present study.

2. Methods

2.1. Single-Difference Residual Reconstruction

In the short baseline carrier phase differential positioning, the satellite orbit error, satellite clock error, receiver clock error, ionospheric error, and tropospheric error can be considered negligible after the difference between stations and satellites. The carrier phase double difference observation equation can be expressed as:

$$\nabla \Delta {\mathsf{\Phi }}_{rb}^{jk}=\nabla \Delta {\rho }_{rb}^{jk}+\lambda \nabla \Delta N_{rb}^{jk}+\nabla \Delta \delta +\nabla \Delta \epsilon$$

(1)

$\nabla \Delta$ denotes the double difference operator, $\rho$ represents the geometric distance between satellite and receiver, $N$ is the carrier phase integer ambiguity, $\delta$ indicates the multipath error, $\epsilon$ refers to the observation noise.

In the static short baseline environment, multipath error has a periodic character. We can employ sidereal filtering to correct the multipath error. The sideral filtering algorithm applying the double-difference residual is more mature. However, it is necessary to solve the problem of reference satellite change or satellite loss in the process of using double difference. Single-difference residual does not have so many restrictions; it only involves the residual of the difference between stations. In addition, when the multipath correction is carried out, it only needs to refer to the repetition period of its own satellite multipath. Therefore, this contribution converts the double-difference residual to the single-difference residual.

Based on obtaining the double-difference residual, the “zero differences” (ZDs) theory proposed by Alber is adopted to obtain the single-difference residual [34]:

$$\mathit{D}\mathit{s}=\mathit{d}$$

(2)

where $\mathit{D}$ is the coefficient matrix, $\mathit{s}$ is the vector composed of single-difference residual, and $\mathit{d}$ is the vector composed of double-difference residual. In order to ensure that the $\mathit{D}$ matrix is invertible, additional constraints need to be attached to the single-difference residual:

$$\left[\begin{array}{ccccc}{\omega }_1& {\omega }_2& {\omega }_3& \cdots & {\omega }_n\\ 1& -1& 0& \cdots & 0\\ 1& 0& -1& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& 0& 0& \cdots & -1\end{array}\right]\left[\begin{array}{c}{s}_{rb}^1\\ s_{rb}^2\\ s_{rb}^3\\ ⋮\\ s_{rb}^n\end{array}\right]=\left[\begin{array}{c}{\omega }_1{s}_{rb}^1+\cdots {\omega }_n{s}_{rb}^n\\ s_{rb}^1-s_{rb}^2\\ s_{rb}^1-s_{rb}^3\\ ⋮\\ s_{rb}^1-s_{rb}^n\end{array}\right]=\left[\begin{array}{c}{\displaystyle \sum {\omega }_i{s}_{rb}^i}\\ d_{rb}^{12}\\ d_{rb}^{13}\\ ⋮\\ d_{rb}^{1n}\end{array}\right]$$

(3)

where ${\omega }_i$ is the weight coefficient of the satellite $i$ on the reference station and the user station, and ${\displaystyle \sum {\omega }_i{s}_{rb}^i}=0$ is set to ensure that the matrix is invertible. It is generally believed that the magnitude of multipath error is affected by the elevation angle of the satellite. Additionally, the weight coefficient ${\omega }_i$ is determined by the elevation angle of satellite $i$, and the weight can be determined by ${\omega }_i={\sin }^2\left({\theta }_i\right)$ where ${\theta }_i\in \left(0,\pi /2\right)$, for $i=1, 2, \dots \mathrm{n}$.

2.2. Tikhonov Regularization Modeling

2.2.1. Tikhonov Regularization Modeling Extraction Multipath

According to the above, the transformed single-difference residual is the original data for multipath modeling, composed of multipath error $\delta$ and observation noise $\epsilon$. It can be expressed as:

$${\varphi }_i=m_i+{\epsilon }_i\hspace{1em}i=1 , 2 , \dots \mathrm{n}$$

(4)

where the subscript $i$ denotes the observation epoch, $\varphi$ represents the transformed single-difference residual sequence, $m$ is the multipath sequence to be extracted, and $\epsilon$ refers to the observation noise. In order to separate multipath sequences, the regularization objective function can be constructed as follows:

$$J\left({\mathrm{m}}_i\right)={\displaystyle \sum _{i=1}^n{\omega }_i{\left({\varphi }_i-m_i\right)}^2}+\alpha {\displaystyle \sum _{i=1}^{n-1}{\left(m_{i+1}-m_i\right)}^2}$$

(5)

where ${\omega }_i$ is the weight coefficient determined by the elevation angle of the satellite, and $\alpha$ is the regularization parameter. Equation (5) is the first-order difference model of Tikhonov regularization. Different forms of objective functions have been discussed in existing literature [25], and there is no obvious difference in their effects. Therefore, only the first-order difference model is utilized in this contribution. Equation (5) can be expressed in matrix form as:

$$\mathit{J}={\left(\mathit{\varphi }-\mathit{m}\right)}^{\mathrm{T}}\mathit{W}\left(\mathit{\varphi }-\mathit{m}\right)+\alpha {\mathit{m}}^{\mathrm{T}}\mathit{R}\mathit{m}$$

(6)

where $\mathit{W}$ denotes the weight matrix, $\mathit{R}$ represents the regularization matrix, and $\alpha$ is the regularization parameter.

To obtain the multipath sequence in Equation (6), its minimizer satisfies the following necessary conditions:

$$\left(\mathit{W}+\alpha \mathit{R}\right) \hat{\mathit{m}}=\mathit{W}\mathit{\varphi }$$

(7)

Note that the coefficient matrix $\mathit{W}+\alpha \mathit{R}$ is invertible only if the null spaces of $\mathit{W}$ and $\mathit{R}$ have a trivial intersection, i.e., $Null\left(\mathit{W}\right)\cap Null\left(\mathit{R}\right)=0$.

2.2.2. Determination of the Regularized Matrix $\mathit{R}$

The regularization matrix is determined by selecting the sum of squares of the difference between the errors of the two adjacent multipath models, that is:

$${\mathit{m}}^{\mathrm{T}}\mathit{R}\mathit{m}={\displaystyle \sum _{i=1}^{n-1}{\left(m_{i+1}-m_i\right)}^2}$$

(8)

the matrix $\mathit{R}$ results

$$\mathit{R}={\mathit{G}}^{\mathrm{T}}\mathit{G}=\left[\begin{array}{rrrrr}\hfill 1& \hfill -1& \hfill 0& \hfill \cdots & \hfill 0\\ \hfill -1& \hfill 2& \hfill -1& \hfill & \hfill \\ \hfill & \hfill \ddots & \hfill \ddots & \hfill \ddots & \hfill \\ \hfill & \hfill & \hfill -1& \hfill 2& \hfill -1\\ \hfill 0& \hfill \cdots & \hfill 0& \hfill -1& \hfill 1\end{array}\right]$$

(9)

where

$$\begin{array}{l}\underset{(\mathrm{n}-1)\times \mathrm{n}}{\mathit{G}}=\left[\begin{array}{ccccc}-1& 1& & & \\ & -1& 1& & \\ & & \ddots & & \\ & & -1& 1& \\ & & & -1& 1\end{array}\right]\\ \end{array}$$

(10)

$\mathit{R}$ is the discretization of the Laplacian with Neumann boundary conditions and has null space, $Null\left(\mathit{R}\right)=\left\{1\right\}$, where $1={\left[\begin{array}{cccc}1& 1& \cdots & 1\end{array}\right]}^T$, see [35]. Therefore, $Null\left(\mathit{W}\right)\cap Null\left(\mathit{R}\right)=0$ ensures the invertibility of the matrix $\mathit{W}+\alpha \mathit{R}$ only if ${\displaystyle \sum {\omega }_i}\ne 0$. This condition is satisfied since the elevation angle is $\left(0,\pi /2\right)$, then ${\omega }_i>0$.

The corresponding multipath sequence can be estimated by Equation (7). In the actual processing process, in order to improve the efficiency of data processing, matrix decomposition algorithms such as the Thomas algorithm can be used to simplify the equation [36,37].

2.2.3. Determination of Regularization Parameter $\alpha$

The regularization parameter is used to adjust the balance between the residual and constraint terms. The selection of the regularization parameter is sensitive to the quality of the solutions. Therefore, it is particularly important to select an appropriate regularization parameter. The bootstrap method evaluated the multipath modeling error to obtain the regularization parameter [38]. The candidate set of regularization parameters is determined in this study, we specify $\alpha =\left\{0.1,1,10,50,100,\dots \right\}$. For each regularization parameter $\alpha$ in the candidate set, its corresponding multipath sequence is determined separately by Equation (7). The following normalized residual (residual of the carrier phase single-difference residual) is defined as:

$${\eta }_i={\omega }_i\left({\varphi }_i-{\widehat{m}}_i\right)\hspace{1em}\mathrm{i}=1 , 2 , \dots \mathrm{n}$$

(11)

where ${\eta }_i$ is the residual of the carrier phase single-difference residual; ${\omega }_i$ denotes the weight determined by the satellite elevation angle; ${\varphi }_i$ is the carrier phase single-difference residual; and ${\widehat{m}}_i$ is the extracted multipath error.

A new sequence $\left\{{\widehat{\eta }}_i\right\}$ is obtained from a random sample with put-back in the normalized residual sequence ${\eta }_i$. Adding $\left\{{\widehat{\eta }}_i/{\omega }_i\right\}$ to $\left\{\widehat{m}\right\}$ generates a resampled carrier phase single-difference residual vector ${\widehat{\varphi }}_i$. Repeat this resampling process $B$ times ($B$ is taken 100 times in this paper), using the carrier phase single-difference residual vector of the $b$th $\left(1\le b\le B\right)$ resampling. The corresponding estimate can be represented as ${\widehat{m}}_b\left(\alpha \right)$. Taking the average as the final estimate of the regularization parameter, namely:

$$\widehat{m}\left(\alpha \right)=\frac{1}{B+1}{\displaystyle \sum _{b=0}^B{\widehat{m}}_b}\left(\alpha \right)$$

(12)

When the subscript $b=0$, $\varphi$ is the original carrier phase single-difference residual. Define the modeling error as follows:

$$E\left(\alpha \right)=\frac{1}{nB}{\displaystyle \sum _{b=0}^B{\left[{\widehat{m}}_b\left(\alpha \right)-\widehat{m}\left(\alpha \right)\right]}^{\mathrm{T}}\left[{\widehat{m}}_b\left(\alpha \right)-\widehat{m}\left(\alpha \right)\right]}$$

(13)

The regularization parameter is determined as

$$\widehat{\alpha }=argmin E\left(\alpha \right)\hspace{1em}\alpha \in \left\{0.1 , 1 , 10 , 50 , 100 , \dots \right\}$$

(14)

In order to further optimize the regularization parameter, this study proposes a method for the circular selection of regularization parameters. Based on the regularization parameter selected by the bootstrap method, the regularization parameter interval $\left[\alpha -m\alpha ,\alpha +n\alpha \right]$ is generated ($m$ and $n$ are positive numbers), and multiple regularization parameters ${\alpha }_i$ are selected cyclically in the interval in steps of 0.1 $\alpha$. The multipath error sequence corresponding to each ${\alpha }_i$ is extracted separately. By exploring the relationship between modeling error and multipath error, the multipath error is extracted to the maximum extent and the modeling error is minimized. The corresponding regularization parameter is taken as the optimal regularization parameter. The specific process is:

• The bootstrap method is used to select the regularization parameter $\alpha$.
• Substitute into Equation (7) and calculate the corresponding multipath error $\widehat{m}$.
• Compute each ${\alpha }_i$ in the interval $\left[\alpha -m\alpha ,\alpha +n\alpha \right]$ cyclically in steps of 0.1 $\alpha$. According to our experience, when m equals 0.1 and n equals 2, the computational efficiency is within an acceptable range while the optimal regularization parameter is guaranteed.
• Repeat step (2) and calculate the multipath error ${\widehat{m}}_i$ and modeling error $E\left({\alpha }_i\right)$ corresponding to each ${\alpha }_i$.
• Output the regularization parameter ${\alpha }_{opt}$ corresponding to $min\left\{E\left({\alpha }_i\right)\right\}$.

2.3. Wavelet Filtering

The wavelet transform analyzes the signal on a time scale and has the characteristics of multi-resolution multi-analysis [39]. Considering the excellent time-domain and frequency-domain characteristics of wavelets, the extraction and elimination of multipath signals using the wavelet method have great superiority. In this study, the effectiveness and superiority of the proposed multipath mitigation strategy are verified by visual comparison with the wavelet filtering method.

Let the basis wavelet be $\psi \left(t\right)$, and the wavelet function generated by the basis wavelet is expressed as:

$${\psi }_{a,b}\left(t\right)={\left|a\right|}^{-\frac{1}{2}}\psi \left(\frac{t-b}{a}\right)a\in R^{+};b\in R$$

(15)

where $a$ and $b$ are the scale and translation parameters, respectively.

Since the carrier phase single-difference residual sequences used are discrete, it is necessary to discretize $a$ and $b$. The corresponding discrete wavelet function of Equation (14) is expressed as:

$${\psi }_{j,k}\left(t\right)=a_0{}^{-\frac{j}{2}}\psi \left(a_0{}^{-\frac{j}{2}}t-k{b}_0\right)\hspace{1em}\mathrm{j},\mathrm{k}\in \mathrm{Z}$$

(16)

where $a_0$ and $b_0$ are the scale parameter and translation parameter of the discrete wavelet function, respectively.

The coefficients of the discrete wavelet transform can be expressed as:

$$C_{j,k}\left(t\right)={\displaystyle {\int }_{-\infty }^{\infty }f\left(t\right)}{\overline{\psi }}_{j,k}\left(t\right)dt$$

(17)

where $f\left(t\right)$ is the sequence of carrier phase single-difference residual to be denoised.

Its reconstruction equation is:

$$f\left(t\right)=c{\displaystyle \sum _{-\infty }^{\infty }{\displaystyle \sum _{-\infty }^{\infty }{W}_{j,k}}}{\psi }_{j,k}\left(t\right)$$

(18)

where $c$ is a constant independent of the signal.

2.4. Attitude Determination

The three attitude angles, namely yaw, pitch, and roll, describe the direction of the carrier coordinate frame relative to the navigation coordinate frame [40]. In this study, we only consider the case of dual antennas and no longer consider the roll angle. The two antennas are assumed to be arranged sequentially on the carrier platform and parallel to the carrier direction. Taking the back antenna as the origin and selecting the direction of the two antennas as the Y-axis, the X-axis, and the Z-axis are perpendicular to the Y-axis and point to the right and the zenith, respectively. The conversion relationship of the baseline vector between the two coordinate frames is:

$$\begin{array}{c}{\left[E\hspace{1em}N\hspace{1em}U\right]}^{\mathrm{T}}=R_b^l·{\left[0\hspace{1em}l\hspace{1em}0\right]}^{\mathrm{T}}\\ R_b^l=\left[\begin{array}{ccc}\cos \gamma \cos \phi +\sin \gamma \sin \theta \sin \phi & \cos \theta \sin \phi & \sin \gamma \cos \phi -\cos \gamma \sin \theta \sin \phi \\ -\cos \gamma \sin \phi +\sin \gamma \sin \theta \cos \phi & \cos \theta \cos \phi & -\sin \gamma \sin \phi -\cos \gamma \sin \theta \cos \phi \\ -sin\gamma \cos \theta & \sin \theta & \cos \gamma \cos \theta \end{array}\right]\end{array}$$

(19)

where $l$ is the a priori baseline length, $R_b^l$ denotes the coordinate transformation matrix, and $\phi$, $\theta$ and $\gamma$ represent the yaw, pitch, and roll angle, respectively.

According to Equation (19), when the yaw and pitch angle is known, the baseline vector can be calculated as:

$${\left[\begin{array}{ccc}E& N& U\end{array}\right]}^{\mathrm{T}}={\left[\begin{array}{ccc}\cos \theta \sin \phi & \cos \theta \cos \phi & \sin \theta \end{array}\right]}^{\mathrm{T}}·l$$

(20)

Similarly, the yaw angle $\phi$ and pitch angle $\theta$ can be calculated using the known baseline vector:

$$\begin{array}{l}\phi =\arctan \left(E/N\right)\\ \theta =\arctan \left(U/\sqrt{E^2+N^2}\right)\end{array}$$

(21)

3. Data Processing and Analysis

In order to verify the consistency of the orbit repetition period and multipath repetition period of the Beidou satellite and the feasibility of the regularization method, this study reconstructs the single-difference residual, extracts the multipath model of each satellite, and evaluates the positioning and attitude results. The experimental data are obtained from the short baseline data from the GNSS Research Center at Curtin University. The observation data collected from CUT0 and CUTC stations are used in this paper to validate the performance of the proposed method. All antennas used are deployed on the roof of the building, and the actual location distribution and antenna information of all stations are presented on the website [41]. The station distance is 7.99 m, the receiver model is TRIMBLE NETR9, the time is from 18 August 2022 (DOY-230) to 31 August 2022 (DOY-243), the sampling interval is 30 s, and the elevation mask is set to 10°. The implementation process of the proposed method is based on the GNSS data processing software RTKLIB and the MATLAB platform [42]. The flow of the multipath error reduction strategy is shown in Figure 1. In static mode, the data from the previous period are obtained for processing. The final baseline vector is estimated after the ambiguity is fixed, and the double-difference residual is estimated for each satellite by substituting into the carrier phase double difference equation epoch by epoch. This approach considers the need to resolve reference satellite changes or losses during the use of double difference. Therefore, we convert the estimated double-difference residual into a single-difference residual, and the reconstructed single-difference residual contains multipath error and other random noise.

The BDS, GEO, and IGSO satellites have good repeatability in their orbits on two adjacent days due to their unique geosynchronous operation. The corresponding orbit repetition period of these two types of satellites is about one sidereal day. The orbits of the MEO satellite have no repetition within two adjacent days. In contrast, there is a strong correlation between the orbit of this type of satellite after an interval of 7 days. Therefore, its corresponding orbital repetition period is about 7 sideral days. Figure 2 shows the carrier phase single-difference residual of partial satellites at DOY-237–238 (DOY-231–238). It can be seen from the figure that for the GEO and IGSO satellites in DOY-237 and DOY-238, the carrier phase single-difference residual has great similarity. This similarity is shown by the fact that the single-difference residual curves of the second day exhibit repeatability with a certain time advance compared to the single-difference residual curves of the first day. In the same way, the MEO satellite has great similarity in the carrier phase single-difference residual in DOY-231 and DOY-238. However, the corresponding time advance is obviously different in order of magnitude from the GEO and IGSO satellite.

As described in the introduction to Section 1, one of the keys to sidereal filtering is the accurate estimation of multipath repetition time. In this study, based on the ORTM, the day is divided into 24 periods according to the broadcast ephemeris, and the MRT is estimated in segments using the orbital parameters of each hour and corrected in segments. The data analysis shows that the MRT estimated using the segmentation idea is slightly better than the ORTM method, consistent with the reference conclusion [30]. Figure 3 shows the mean MRT estimates for all satellites and the MRT estimates for each segment. The yellow bars indicate the corresponding mean MRT estimates for the GEO and IGSO satellites, and the purple bars indicate the corresponding mean MRT estimates for the MEO satellite. The red color indicates the MRT estimates corresponding to each segment for GEO and IGSO satellites, and the green color indicates the MRT estimates corresponding to each segment for the MEO satellite. The MRT of the single-difference residual of the GEO and IGSO satellites is different in the corresponding period, which differs only by tens of seconds, with an average MRT of about 239 s. For the MEO satellite, the average MRT corresponding to the single-difference residual is about 1698 s in the corresponding period.

In the regularized extraction process, the squared residual term is not robust to outliers due to its inclusion in the objective function. Therefore, cubic spline interpolation was used to pre-interpolate the single-difference residual to replace the outliers and fill in the missing data [43]. High-quality data are provided for subsequent extraction of multipath. According to the Pauta criterion (3σ criterion), if the deviation between the data value and the interpolation is three times the standard deviation of the corresponding carrier phase measurement, it is regarded as an outlier. Additionally, the data not exceeding the limit remain unchanged [44]. We denote the Tikhonov regularization method using the bootstrap method to select the regularization parameter as the TB method and the Tikhonov regularization method using a circular selection of the regularization parameter as the TC method. Figure 4 displays the multipath error extracted by C01 (left) and C07 (middle) at DOY-237 and C28 (right) at DOY-231 using TB (top) and TC (bottom) methods. It can be seen from Figure 4 that it is feasible to extract multipath error correctly by Tikhonov regularization, and multipath error is more smooth than the original carrier phase single-difference residual. The multipath error extracted using the TC method is slightly better than the TB method. This means that the multipath error extracted by obtaining the optimal regularization parameter as much as possible is more accurate when the multipath error is extracted using the Tikhonov regularization method.

To reflect the validity of the method, the wavelet soft threshold denoising method is also used for data analysis to visually compare the relative performance of the proposed method. The selected wavelet basis function is sym6, the number of decomposition layers is 4, and the filtering threshold uses the default threshold provided by MATLAB. The selection of the wavelet basis function and the number of decomposition layers can be found in the literature [45,46,47,48]. Figure 5 displays the multipath error extracted by the C01 (left) and C07 (middle) at DOY-237 and the C28 (right) at DOY-231 using the TC method (top) and the wavelet filtering method (bottom). From Figure 5, it can be concluded that the proposed method has a good denoising effect. Additionally, the systematic multipath error is extracted, which is consistent with the residual volatility and has good fidelity. Meanwhile, compared with classical wavelet filtering, the proposed method extracts multipath error curves that are similar to wavelet filtering, and the improvement of phase residual before and after filtering is better than wavelet filtering.

After the above process, the two key steps of sideral filtering are completed: the calculation of MRT and the establishment of a multipath error extraction model. Therefore, it is possible to apply the previous period’s data to the current period’s data. The multipath error established in the previous period is subtracted from the current period using the calculated time advance. The carrier phase single-difference residual before and after applying the sideral filtering of the TB and TC methods to C04, C06, and C19 at DOY-238 are shown in Figure 6. It can be observed that the multipath error has been effectively suppressed with significant effect. The filtered a posteriori single-difference residual (cyan line) has apparently removed the significant trend in the original single-difference residual and exhibited random characteristics, with the better part of the effect approaching white noise. Figure 7 displays the carrier phase single-difference residual before and after the sideral filtering of the TC method and the wavelet filtering method applied to C04, C06, and C19 at DOY-238. It can be observed from Figure 7 that the regularization method is slightly better than the wavelet method, which further validates the effectiveness of the proposed method.

The RMS and the amount of improvement in the carrier phase single-difference residual before and after multipath error mitigation by sidereal filtering are shown in Figure 8. We verify the superiority of the proposed TC method by comparing the average percentage of RMS improvement in the single-difference residual sequence before and after sidereal filtering for all visible satellites. This is defined as the average RMS improvement rate. It can be observed from the figure that the average RMS improvement rate for applying TB, TC, and wavelet filtering methods to all visible satellites is 33.9%, 40.5%, and 32.6%, respectively. In response to the unique constellation characteristics of BDS, we counted the improvement of each type of satellite. For the GEO satellite, adopting TB, TC, and wavelet filtering methods, the average RMS improvement rate is 37.9%, 45.9%, and 36.7%, respectively. For the IGSO satellite, the average RMS improvement rate of the three methods is 30.5%, 38.2%, and 26.4%, respectively. For the MEO satellite, the average RMS improvement rate is 33.3%, 37.5%, and 34.5% after applying TB, TC, and wavelet methods, respectively. It can be concluded that the proposed method has the best improvement effect after applying sidereal filtering. In addition, the multipath improvement for the GEO satellite is more significant than for the other two types of satellites.

The multipath error has a significant effect on the success rate of ambiguity fixing. The improvement in the success rate of ambiguity fixing after multipath mitigation will further prove the superiority of the proposed method. The success rate of ambiguity fixing and RMS before and after multipath error improvement are shown in

Table 1. The success rate of ambiguity fixing and RMS statistics before and after multipath error correction (unit: mm).

Type	Uncorrected	After Correcting the Multipath
		TB Method	Imp.	TC Method	Imp.	Wavelet Filtering	Imp.
Fixed Rate	89.4%	97.3%	7.9%	97.6%	8.2%	97.5%	8.1%
E	16.5	12.7	23%	12.4	24.8%	12.8	22.4%
N	3.8	2.8	26.3%	2.8	26.3%	2.9	23.7%
U	9.6	5.6	41.7%	5.5	42.7%	5.6	41.7%

. After sidereal filtering multipath mitigation, the success rate of ambiguity fixing for TB, TC, and wavelet filtering methods improved by 7.9%/8.2%/8.1%. The positioning accuracy is higher in the N and U components than in the E component, and the positioning results are satisfactory. The correction effect of the three methods is very significant. After the TB, TC, and wavelet filtering methods, the E/N/U component positioning accuracy improved to 23%/26.3%/41.7%, 24.8%/26.3%/42.7%, and 22.4%/23.7%/41.7%, respectively. In general, it is more effective to correct multipath errors using the TC method.

Table 2. Coordinate representation of the two antennas in the ECEF coordinate frame (unit: m).

Antennas	X	Y	Z
CUT0	−2364337.6799	4870285.6506	−3360809.3985
CUTC	−2364331.4902	4870284.8979	−3360814.3954

lists the known coordinates of the two antennas in the ECEF coordinate frame. After further coordinate conversion, the coordinates of the baseline vector in the carrier coordinate frame are easily obtained as (−5.2339, −6.0303, −0.2133), and the baseline length is 7.99 m. Based on this, it is easy to calculate the reference attitude information for this experiment: the yaw and pitch angles are 40.9558° and −1.5302°, respectively.

Table 3. Statistics of the average value of attitude angle of TC method in each hour (unit: deg).

Attitude	0–1 h	1–2 h	2–3 h	3–4 h
Yaw	40.9972	40.9659	40.9627	40.9574
Pitch	−1.4953	−1.5345	−1.5272	−1.5351

shows the average value of the attitude angle for each hour after extracting the multipath by the TC method. It can be observed that the average value of the attitude angle in the first hour differs greatly from the reference value. With the increased epoch, the calculated attitude angle gradually approaches the reference value. The possible reason for the above results is that forward filtering is used and the ambiguity of some epochs is not fixed at the beginning. As the number of epochs increases, the number of fixed ambiguities increases. The attitude resolution accuracy will remain relatively stable.

The attitude angle error variation using the original data and after applying the TC method is shown in Figure 9. It is easily observed from Figure 9 that the difference between the attitude and the reference value before and after the sidereal day filtering mitigation multipath, i.e., the attitude angle error, is concentrated at 0°. The change in attitude angle is not very obvious from the error plot of the time series. In order to be able to show more clearly the changes in the attitude angle before and after mitigating the multipath using the TC method. We use the cumulative distribution function of the absolute value of the error in Figure 10 to further illustrate. Through further data analysis, the RMS of the yaw angle is within 0.03° and the RMS of the pitch angle is about 0.01° after applying the TC method to weaken the multipath. We evaluate the effect on attitude angle after adopting the TC method to mitigate multipath and the RMS improvement rate of attitude angle close to the reference attitude angle before and after multipath mitigation. Compared with the attitude angles resolved from the original observations, the RMS improvement rate of yaw and pitch angles after mitigating the multipath using the TC method is 22.9% and 12.6%, respectively.

4. Discussion

Sidereal filtering is a common method to deal with multipath error, which is basically based on the fact that the periodic motion of satellites causes multipath error to show periodic variations. Currently, the main methods to mitigate multipath error using periodic repeatability include correction based on the observation domain and correction based on the coordinate domain. It has been confirmed in the literature using GPS data that the results based on the observation domain are significantly better than the results of multipath error correction based on the coordinate domain [24]. When applying the multipath error correction in the coordinate domain, only the mean orbit repetition period of the observed satellites is considered, ignoring the differences between satellites. This method only applies to GPS and GLONASS with the same orbit but cannot be effectively applied to the BDS multipath effect with three different types of orbits. Therefore, this study corrects the multipath error based on the observation domain. It obtains the double-difference residual sequence through the static resolution of the previous period. Subsequently, it performs the regularization and wavelet filtering methods to denoise the converted single-difference residual. The multipath correction model of each satellite is then derived to correct the carrier phase residual of the current period. Finally, this study evaluates the performance of the proposed method using the positioning and attitude results. The method is highly adaptable and achieves good results by considering the repetition period of each satellite through the segmentation idea, which is complemented by efficient numerical calculation methods such as the Thoms algorithm. This makes the computation within an acceptable range.

In the experiment, it can be seen in Figure 2 that the multipath error is systematic. This means that multipath errors exist, and their effects can be weakened by building the corresponding mathematical models. In addition, the strong similarity of adjacent periods suggests the feasibility of modeling the previous period’s data to correct the current period’s data. Figure 3 quantitatively depicts the repetition period of each satellite of the BDS. It can be seen that the MRT corresponding to each hour within the corresponding period varies by tens of seconds. Therefore, the MRT estimated using the segmented orbital parameters per hour is more accurate. In Figure 4 and Figure 5, the multipath extraction modeling using regularization and wavelet filtering methods shows good performance. This confirms the effectiveness of multipath extraction modeling with the regularization method. In Figure 6, Figure 7, Figure 8 and Figure 9, after applying the multipath extraction modeling established in this study to the sidereal filtering, the multipath error is obviously suppressed. The improvement effect of the proposed method on multipath is better than that of the wavelet filtering method. This can be further verified in the statistical results of positioning and attitude measurement in Table 1 and Table 3.

In summary, the Tikhonov regularization method based on parameter optimization proposed in this study has high modeling accuracy and numerical computation efficiency and provides an optional method for GNSS multipath extraction and reduction.

5. Conclusions

In static or quasistatic modes, the multipath error is to some extent systematic rather than purely random and should be considered a signal rather than noise. Based on this fact, the multipath error can be empirically modeled using the repeatability feature. This contribution confronts the unique operation of BDS based on the reconstructed carrier phase single-difference residual and extracts the multipath error sequence using parameter-optimized Tikhonov regularization to correct the multipath error for each BDS satellite separately. Additionally, the correctness and effectiveness of the method are verified by experiments.

Accurate estimation of multipath repetition time is one of the key steps in sidereal filtering. The MRT estimated using the hourly orbit parameter segments in the broadcast ephemeris is slightly better than the ORTM method. The intuitive performance comparison between the regularization method and the wavelet soft thresholding denoising method shows that it is feasible to correctly extract the multipath error by adopting Tikhonov regularization. Furthermore, the multipath error is smoother compared to the original residual. Moreover, the multipath error extracted employing the TC method proposed in this study is better than the TB method, which implies that the optimal regularization parameter should be estimated as much as possible when applying the Tikhonov regularization method. For the carrier phase single-difference residual, the average RMS improvements of GEO, IGSO, and MEO satellites are 45.9%, 38.2%, and 37.5%, respectively. Using the proposed method in this study, the multipath improvement effect for GEO satellites is more significant. For coordinate residual, the positioning accuracy is improved by 24.8% in the east–west direction, 26.3% in the north–south direction, and 42.7% in the elevation direction after multipath error mitigation. After further attitude angle resolution, the yaw angle accuracy is improved by 22.9%, and the pitch angle accuracy is improved by 12.6%. As a final note, it is necessary to point out that although this study uses parametric optimization and Tikhonov regularization to mitigate the multipath error of BDS, it is also possible to incorporate a more accurate regularization matrix R, as proposed in [35]. This approach may provide better results.