## Abstract

To analyze the unmodeled systematic errors in BDS, a filter-assisted Partly Ensemble Empirical Mode Decomposition (PEEMD) method is combined with Hilbert spectrum analysis to extract the feature information from BDS preprocessed orbit determination residuals (PODR). The results show that the feature extraction method can effectively extract a period of about 1 day for GEO/IGSO satellites and a period of about 13 h for MEO satellites. The results of the chi-square test show that the remaining PODR follow a normal distribution. Statistical results from the 3-day experiment indicate that the application of the extracted feature information for BeiDou-only SPP improves positioning accuracy by 20%, 19%, and 23% in the east, north, and upward directions, respectively. BeiDou-only PPP experiments show that the application of the extracted feature information reduces the static PPP convergence time by 34%, 10%, and 21% and the kinematic PPP convergence time by 25%, 11%, and 9% in three coordinate components, respectively.

## 1 INTRODUCTION

Since the end of 2012, the BeiDou Navigation Satellite (BDS) system has realized regional positioning services over the Asia-Pacific region. Its current heterogeneous constellation consists of five Geostationary Orbit (GEO) satellites, six Inclined Geosynchronous Orbit (IGSO) satellites, and three Medium-altitude Earth Orbit (MEO) satellites. The distribution of monitoring stations and the special constellation of BDS mean that its orbit determination is inferior to that of the global positioning system (GPS). In particular, the GEO satellites suffer from near-static viewing geometry with respect to the receivers, resulting in the worst orbit accuracy among the BDS satellite.^{1}

The orbit determination residuals (ODR) are a manifestation of precision, and the orbit determination accuracy can be evaluated by analyzing the ODR.^{2} Improving the accuracy of the satellite orbit, clock offset, and relevant model parameters is an effective method of enhancing BDS positioning, navigation, and timing performance.^{3} However, unmodeled systematic errors such as multipath effects,^{4,5} ionospheric high-order errors, remaining tropospheric delay, and satellite orbit significantly affect the BDS orbit determination. The extraction of feature information from BDS ODR can help to separate the unmodeled errors for BDS orbit determination. Moreover, the feature information is beneficial in analyzing the causes of the unmodeled systematic errors and improving the precision of BDS satellite orbits.

In order to get orbit determination residuals of several whole orbit periods, the original ODR must be preprocessed, called the preprocessed ODR (PODR). BDS PODR are typically nonlinear and nonstationary and can be studied using time-frequency analysis. Empirical Mode Decomposition (EMD) is an adaptive time-series analysis method.^{6} A nonstationary and nonlinear signal can be decomposed into several Intrinsic Mode Function (IMF) components using EMD. However, one of the primary problems of EMD is the phenomenon of mode mixing caused by intermittence or noise. To resolve this problem, Wu and Huang^{7} proposed a noise-assisted data analysis method called Ensemble EMD (EEMD) by adding white noise to the initial data. Unfortunately, EEMD also has some disadvantages. For example, the IMFs derived from EEMD will inevitably be polluted by the added noise. To effectively eliminate residual noise from the IMFs and restrain the reconstruction residue, Yeh et al^{8} proposed complementary EEMD (CEEMD) to improve the efficiency of the original noise-assisted method by adding noise in positive and negative pairs. However, both EEMD and CEEMD are time-consuming when the number of trials large. Moreover, the IMFs derived from EEMD and CEEMD often include some false and meaningless components. To overcome the disadvantages of EEMD and CEEMD, the noise-assisted Partly Ensemble EMD (PEEMD) approach was proposed.^{9} PEEMD decomposes signals as in CEEMD for the first *j* IMFs and detects the intermittency or noise signal in the original using permutation entropy. Using EMD directly, PEEMD then decomposes the remaining signal. PEEMD avoids the unnecessary ensemble and averaging of all IMFs. Hence, PEEMD significantly reduces the calculation time and improves the accuracy of the resulting IMFs. However, determining the two key parameters, that is, the amplitude of added white noise and the number of ensemble members, is not clearly defined in PEEMD.

On the basis of the original PEEMD, we present a filter-assisted PEEMD method in which the two key parameters can be determined by filtering the original signals. Then, combining the filter-assisted PEEMD with Hilbert spectrum analysis, a feature extraction method is proposed for BDS PODR. The intermediate- and low-frequency feature information is then used to establish the model signal. The remaining PODR, which are the differences between the original signal and the model signal, are also established. To validate the effectiveness of the feature extraction method, the distribution characteristics of the remaining PODR are quantitatively analyzed using a chi-square test. Furthermore, the model signal is applied to BeiDou-only single point positioning (SPP) and static and kinematic precise point positioning (PPP) to further validate the effectiveness of the proposed method.

## 2 METHODOLOGY

BDS PODR are a time series of the differences between the observations and theoretically calculated values, and these directly reflect the precision of orbit determination. PEEMD and Hilbert spectrum analysis are efficient time-series analysis methods, especially for nonstationary and nonlinear time series. Thus, BDS PODR can be analyzed with PEEMD and Hilbert spectrum analysis. The following provides an introduction to feature extraction based on PEEMD and Hilbert spectrum analysis.

### 2.1 PEEMD and Hilbert spectrum analysis

Similar to CEEMD, PEEMD has the advantage of reducing mode mixing and decreasing the reconstruction errors. Moreover, PEEMD can greatly reduce the computational load of CEEMD and improve the accuracy of the resulting IMFs. PEEMD is described as follows^{9}:

1. Suppose

*j*= 1 and add the white noise series {*n*_{i}(*t*)} and {−*n*_{i}(*t*)} to the targeted signal*S*(*t*) as1 2

where *a* represents the amplitude of the added white noise, *i* = 1,2, …, *N*_{e}; *N*_{e} represents the number of pairs of added white noise.

2. Decompose the two signal series and using EMD for the

*j*th IMF mode. Two IMF sets and are obtained, as well as two residue sets and , respectively.Namely,

3 4

3. By constructing an ensemble, the final IMF for the

*j*th rank can be derived by5

4. Calculate the Permutation Entropy (PE)

*P*_{j}of*I*_{j}(*t*). If*P*_{j}is larger than the threshold of PE*θ*_{0}, then set*j*=*j*+ 1 and go to steps (2)–(3); if*P*_{j}is smaller than*θ*_{0}, then go to step (5).5. Separate the first

*j*− 1 IMFs from the original signal and obtain the residue*r*(*t*) by6

6. Decompose

*r*(*t*) completely using EMD:7

7. {

*c*_{k}(*t*)} are regarded as the IMFs following the first*j*− 1 IMFs. The initial signal is given by8

Instead of implementing an ensemble and averaging all IMFs, PEEMD decomposes the original signal in the same way as CEEMD for the first (*j* − 1) IMFs and then uses EMD for the remaining ranks.

After applying a Hilbert transform to each IMF, the signal’s Hilbert energy spectrum can be derived. The Hilbert envelope spectrum is the integral of the Hilbert energy spectrum with respect to time.

The Hilbert transform is defined as

9

where *P* indicates the Cauchy principal value of the integral, *C*_{j}(*t*) denotes the *j*th IMF, and *τ* is the delay time. With the Hilbert transform , the envelope *a*_{j}(*t*), phase *ϑ*_{j}(*t*), and frequency *ω*_{j}(*t*) of the IMF *C*_{j}(*t*) are calculated by^{10}

10 11 12

The time-frequency distribution of the squared amplitude is designated as the Hilbert energy spectrum *H*(*ω, t*).

13

Integrating the Hilbert energy spectrum with respect to time gives the Hilbert envelope spectrum *h*(*f*):

14

The Hilbert envelope spectrum *h*(*f*) offers a measure of the total amplitude (or energy) contribution from each frequency value, providing direct information about the frequency variation of the energy. The input signal can be decomposed into several IMF components by PPEMD, and then the Hilbert transform can be performed for each IMF component to obtain feature frequencies.

### 2.2 Filter-assisted PEEMD

Two important parameters in PEEMD are the amplitude of the added white noise *ε*_{n} and the number of ensemble members *N*_{e}. It is difficult to effectively determine these parameters. The effective high-frequency components are unaffected by the added white Gaussian noise, and this noise changes the distribution characteristics of the extreme intervals of low-frequency components in the signal.^{11} Chen et al^{11} proposed the criterion of adding white Gaussian noise in EEMD^{7}:

15

where *α* represents the ratio of the standard deviation of the added white noise amplitude, *ε*_{n}, to the standard deviation of the original signal amplitude, *ε*_{0}, that is, ; *σ* represents the ratio of the standard deviation of the amplitude of effective components in the original signal, *ε*_{h}, to *ε*_{0}, that is, . Chen et al^{11} suggested that, in most cases, when , the phenomenon of mode mixing can be effectively avoided in the signal decomposition.

The number of ensemble members *N*_{e} is determined by the following formula^{12}:

16

where *e* is the final standard deviation of the error, defined as the difference between the input signal and the corresponding IMFs.

With reference to the above method of determining the two EEMD parameters, a filter-assisted PEEMD method is proposed. The two parameters can be determined by filtering the original signal in the filter-assisted PEEMD method, which can be described as follows:

1. Determine

*ε*_{h}by filtering the original signal using a high-pass or band-pass filter.2. Determine

*ε*_{0}.3. Determine

*σ*and*α*.4. Determine

*N*_{e}.5. Conduct the original PEEMD based on

*α*and*N*_{e}.

The bandwidth of the band-pass filter can be determined according to the characteristics of the original signal; for example, the bandwidth can be set to [3*f*_{s}/8, *f*_{s}/2], where *f*_{s} denotes the sampling frequency. A flowchart of the filter-assisted PEEMD method is shown in Figure 1.

### 2.3 Improved feature extraction method

Based on the combination of the filter-assisted PEEMD method and Hilbert spectrum analysis, a feature extraction method for BDS PODR is proposed. A flowchart of the feature extraction method is shown in Figure 2.

The filter-assisted PEEMD method decomposes the time series of BDS PODR into a finite number of IMFs, *C*_{j}(*t*), *j* = 1, …, *n*, and a residue, *r*_{n}(*t*). Each IMF is the energy associated with various intrinsic time scales extracted from the original signal. The application of the Hilbert transform to each IMF extracts the feature frequencies of BDS PODR. The resulting feature information can be further analyzed in the frequency domain. Based on the feature frequencies, the model signals and remaining PODR can be established.

The chi-square test is an effective method of verifying whether the remaining PODR obey the normal distribution. When the remaining PODR exhibit normal distribution characteristics, it indicates that all feature information has been included in the model signal. BeiDou-only SPP and static and kinematic PPP experiments are carried out to evaluate the effectiveness of the extracted model signals.

## 3 DATA ACQUISITION AND EXPERIMENTS

Determining satellite orbit requires the observations of several stations. Therefore, after the orbit is determined, we can get the original ODR of carrier phase and pseudorange at each station. To get the residuals through several stations, called PODR, the original ODR have to be preprocessed. The proposed method is used to extract the feature information of PODR. Finally, a chi-square test is conducted to verify the effectiveness of the method, and the effectiveness of the extracted feature information is verified by BeiDou-only SPP and PPP experiments.

### 3.1 Data acquisition and preprocessing

The test datasets used for satellite orbit determination were obtained from 53 Multi-GNSS Experiment (MGEX)^{13} stations of the International GNSS Service (IGS)^{14} and 22 Crustal Movement Observation Network of China (CMONOC) stations from July 12 to 15 in 2016. The distribution of stations is shown in Figure 3.

All stations are capable of collecting GPS and BDS observations. The collected datasets were processed by the Position and Navigation Data Analyst (PANDA) software.^{15} Orbit determination method is “one-step” integrated precise orbit determination.^{16} The data-processing observation model and the dynamic model selection are shown in Table 1.

For GEO satellites, we choose a station. For IGSO satellites, we choose two stations, and these two stations are combined for global coverage. For the MEO satellites, the number of stations depend on the stations’ distribution and the continuity of the stations data. For IGSO and MEO satellite, when multiple stations contain the same satellites, the residuals of the observations on the station with higher satellite elevation angle are selected.

After being processed by PANDA, the positions of each satellite are obtained. Simultaneously, the original ODR of carrier phase and pseudorange can be determined at a sampling interval of 300 s. As BDS includes three satellite types, different preprocessing strategies are needed to obtain PODR. The GEO satellites require one station, IGSO satellites require a minimum of two stations, and MEO satellites require more stations. The number of stations depend on the stations’ distribution and the continuity of the stations’ data. Table 2 lists the selected stations for GEO-C003, IGSO-C005, and MEO-C012, respectively. Figure 4 shows the ground tracks of C003, C005, C012, and the locations of the selected stations.

We use 3 days of PODR for analysis. Figure 5 shows the PODR of C003, C005, and C012. It can be seen that the pseudorange PODR are significantly larger than the carrier phase PODR, and all residuals contain a lot of random noise. Although the C012 PODR present disordered features, its modulation information can be found by careful analysis. Additionally, there are several obvious high and low peaks in the PODR of C003 and C005. The BDS PODR do not present random distribution features and include some feature information. If this feature information could be extracted, the unmodeled systematic errors would be separated from the BDS PODR, which are expected to improve the BDS orbit determination precision.

### 3.2 Feature extraction of BDS orbit determination residuals

The feature extraction method is used to analyze the carrier phase and pseudorange PODR of satellites C003, C005, and C012 to extract their feature information. The finite impulse response (FIR) filter is selected, and the bandwidth of the band-pass filter is set to [3*f*_{s}/8, *f*_{s}/2], where *f*_{s} is the sampling frequency of BDS POSR. The desired standard deviation of the signal decomposition error *e* is set to 0.002. According to the results of Zheng et al,^{9} the recommended value of *θ*_{0} = 0.5 is used. With the above predetermined values, the key parameters of the filter-assisted PEEMD for each satellite are determined (see Table 3). The results of the filter-assisted PEEMD of the three satellites are shown in Figures 6–8, where the “signal” represents the original signal of BDS PODR; C1, C2, C3, … denotes IMF components (from high frequency to low frequency); and R represents the trend.

The results for C003 in Figure 6 indicate that a component with a period of about 1 day exists in both the pseudorange and carrier phase PODR (C7). Moreover, a component with a period of about 8 h (1/3 of a day) can also be clearly found (C6). The amplitudes of these components in the pseudorange residuals are 25 cm and 24 cm, respectively, and the amplitudes of the carrier phase residuals are 3 mm and 8 mm, respectively. These two components exhibit similar characteristics.

Figure 7 illustrates that a component with a period of about 1 day exists in the pseudorange PODR (C8) and carrier phase PODR (C7) for satellite C005, respectively. The amplitudes of the 1-day period components are 40 cm and 4 mm, respectively. A component with an 8-h period in the pseudorange residuals can also be found (C7), and its amplitude is 13 cm. In the carrier phase residuals, there is a component with a period of about 12 h (C6), and its amplitude is 3 mm.

From the results for satellite C012 in Figure 8, components with a period of about 1 day (C8) and about 12 h (C7) can be clearly found in both the pseudorange and carrier phase PODR. The amplitudes of the pseudorange residuals are approximately 25 cm for C7 and C8 components, whereas those of the carrier phase residuals are 3 cm and 2 cm, respectively.

In conclusion, the periodic oscillatory components over different time scales can be decomposed from BDS PODR by the filter-assisted PEEMD. To extract the feature information accurately, the IMF components with obvious feature are analyzed by the Hilbert spectrum analysis. The resulting Hilbert envelope spectra for the filter-assisted PEEMD results of C003 are shown in Figure 9.

Figure 9 shows that several different feature frequencies exist in the carrier phase and pseudorange PODR. Moreover, three similar feature frequencies are revealed in both residuals. From the results in Figure 9, a 1-day period does exist in the PODR of GEO satellite C003. Also, the orbital period of GEO satellites is about 1 day. Therefore, it is inferred that BDS PODR include some feature information that varies with the BDS orbital period.

To further verify the above inference, the PODR of the other BDS satellites were analyzed by the proposed feature extraction method. The two most obvious feature frequencies of each satellite are listed in Figure 10.

It can be seen that the two most obvious oscillation periods of GEO satellites’ PODR are about 24 h and 12 h or 24 h and 8 h; those of IGSO satellites are about 24 h and 12 h or 24 h and 8 h, and those of MEO satellites are about 26 h and 13 h.

From the above results, we can find that a period of about 1 day exists in the extracted feature information of all GEO/IGSO satellites, which is in good agreement with their orbital periods. Moreover, a 13-h period exists in the feature information extracted from all MEO satellites, which is also in good agreement with their 773-min revolution period.^{17} Thus, it is reasonable to consider that the extracted feature periods are likely to be associated with BDS orbit periods. In conclusion, the feature information hidden in BDS PODR can be accurately extracted using the proposed feature extraction method.

### 3.3 Effectiveness of the feature extraction method

The original PODR can be divided into the model signal and the remaining PODR. Applying the chi-square test to the remaining PODR can verify whether all feature information has been extracted. The extracted model signal can be applied to BeiDou-only SPP and PPP to further evaluate the effectiveness of the proposed method. First, the feature information of every BDS satellite is extracted. Then, the model signals of all BDS satellites are established with the extracted feature information. Finally, BeiDou-only SPP and static and kinematic PPP experiments are conducted to assess the performance by correcting the raw measurements with the extracted model signals. The test datasets collected at three stations, that is, DAE2, KARR, and LHAZ, on three consecutive days, that is, July12–14, 2016, are used for SPP and PPP experiments.

#### 3.3.1 Chi-square test

The chi-square test can be used to verify the distribution characteristics of the PODR. The statistical degrees of freedom and significance level were set to 6 and 0.05, respectively. The critical value corresponding to this degree of freedom and significance level is 12.5916. Taking the pseudorange PODR of C005 as an example, the model signal, the remaining PODR, and the result of the chi-square test are shown in Figure 11.

The chi-square result in the upper panel of Figure 11 is significantly higher than the critical value. Therefore, the original PODR do not obey a normal distribution. This means that some feature information exists in these residuals. However, the chi-square value of the remaining PODR (lower panel of Figure 11) is significantly lower than the critical value, indicating that the remaining PODR obey a normal distribution and all feature information has been extracted in the model signal.

#### 3.3.2 BeiDou-only SPP

The orbit used in the SPP experiment is a precision orbit, and the observation is only pseudorange. Figure 12 shows the positioning errors at DAE2, KARR, and LHAZ in the horizontal and vertical directions. Figure 13 shows the distribution of the SPP positioning errors in three coordinate components. Table 4 shows the statistical results of the SPP positioning errors.

As can be seen from Figure 12, the uncorrected SPP is distributed within plus or minus 3 m, and the corrected SPP is within plus or minus 2 m in the horizontal pane, slightly better than the uncorrected. In the vertical direction, the corrected SPP is less volatile than the uncorrected SPP. It can be seen from the error distribution of Figure 13 that the positioning errors in the corrected SPP result are more concentrated in the small errors over the uncorrected in coordinate components. Table 4 provides the mean, standard deviation (STD), and root mean square (RMS) of SPP positioning errors in three coordinate components. It can be seen from the table that the corrected result is better than the uncorrected result. The RMS of the corrected SPP is improved 20%, 19%, and 23% over the uncorrected SPP in three coordinate components, respectively.

#### 3.3.3 BeiDou-only static PPP

BeiDou-only static PPP experiments were conducted at the same stations. The positioning errors at DAE2, KARR, and LHAZ are shown in Figure 14. Figure 15 shows the measurement residuals at DAE2, KARR, and LHAZ on July 12, 2016, and Figure 16 gives the RMS statistics of static PPP pseudorange residuals. Table 5 shows the statistical results of convergence time for BeiDou-only static PPP.

The positioning results are considered to have converged when the absolute errors are less than 0.05 m in the east and north directions and less than 0.10 m in the up direction in static PPP. As shown in Figure 14, after being corrected by the model signals, the convergence time of static PPP improves, especially in the east and north directions, with respect to the uncorrected results. Table 5 presents the mean, STD, and RMS of the convergence time of static PPP in three coordinate components. The corrected static PPP produces reductions of 34%, 10%, and 21% in three coordinate components, respectively.

Figure 15 illustrates the residuals of corrected and uncorrected static PPP DAE2, KARR, and LHAZ. It is clear that some feature information is significantly eliminated. The residuals of corrected static PPP show better random noise characteristics with respect to the residuals of uncorrected static PPP. Figure 16 illustrates the RMS error of the residuals for corrected and uncorrected static PPP. This shows that the RMS errors of the residuals of corrected static PPP are lower, especially for GEO satellites.

#### 3.3.4 BeiDou-only kinematic PPP

We carried out a simulated kinematic PPP experiment with previous static experimental data. Figure 17 shows the positioning errors of kinematic PPP in three coordinate components. Table 6 provides the convergence time of kinematic PPP. Figure 18 shows the residuals of C016, C005, and C012 for kinematic PPP. Figure 19 shows the statistical results of each satellite residual for kinematic PPP.

In the kinematic PPP experiment, we define that the positioning results have converged when the absolute positioning errors are less than 0.1 m in the east and north directions and less than 0.2 m in the up direction. As shown in Figure 17, the corrected kinematic PPP results converge more quickly than the uncorrected. Table 6 shows the statistical results of the convergence time of kinematic PPP. It can be seen that the corrected kinematic PPP convergence time is reduced by 25%, 11%, and 9% over the uncorrected in three coordinate components, respectively.

Figure 18 depicts the 3-day pseudorange residuals of kinematic PPP for C016, C005, and C012 at DAE2, KARR, and LHAZ. It can be seen that there are some oscillations in the uncorrected residuals of the kinematic PPP, and the corrected residuals show better random noise characteristics. As can be seen from Figure 19, the RMS of the corrected kinematic PPP residuals is smaller than that of the uncorrected, especially for the GEO satellite.

In conclusion, the proposed feature extraction method is capable of separating the feature information from BDS PODR. Using the feature information to correct the pseudorange observations in BeiDou-only SPP can improve the accuracy. Good convergence performance can be achieved with the application of extracted feature information in BDS-only static and kinematic PPP. In addition, the application of feature information can eliminate some oscillation information that existed in residuals.

## 4 CONCLUSION

An improved feature extraction method, based on the combination of the filter-assisted PEEMD and Hilbert spectrum analysis, has been proposed for BDS PODR. BDS PODR were analyzed with the proposed feature extraction method. A period of about 1 day exists in the PODR of all GEO/IGSO satellites, and a period of approximately 13 h exists for all MEO satellites. These periods are in good agreement with the BDS orbital periods. The results of the chi-square test show that the original signal includes some feature information that can be extracted using the proposed approach. Using the extracted model signal to correct pseudorange observations improves the BeiDou-only SPP accuracy by 20%, 19%, and 23% in east, north, and up directions, respectively. The results of BeiDou-only PPP experiments illustrate that the extracted model signal can be used to improve the static PPP convergence time by 34%, 10%, and 21% and the kinematic PPP by 25%, 11%, and 9% in three coordinate components, respectively. These tests together indicate the effectiveness of the feature extraction method. In conclusion, the feature extraction method can be used to separate the unmodeled errors in BDS orbit determination and improve the BeiDou-only SPP accuracy and the convergence time of BeiDou-only static and kinematic PPP.

Future work should study the causes of different feature periods in BDS PODR at the physical level. Moreover, besides BDS, the unmodeled errors of GLONASS and GALILEO will be examined with case-study experiments.

## HOW TO CITE THIS ARTICLE

Hu G, Ye S, Chen D, et al. Extracting unmodeled systematic errors from BDS orbit determination residuals and application in SPP/PPP. *NAVIGATION-US*. 2020;67:275–289. https://doi.org/10.1002/navi.365

## ACKNOWLEDGEMENTS

We are grateful to the anonymous reviewers for their helpful, constructive suggestions, and comments that significantly improved the paper quality. This work is partially supported by the National Science Fund for Distinguished Young Scholars (no. 41525014), National Natural Science Foundation of China (nos. 41074008, 41431069, and 41604028), Research Fund for the Doctoral Program of Higher Education of China (no. 20120141110025), and Anhui Natural Science Foundation (no. 1708085QD83).

## Footnotes

**Funding information**Anhui Natural Science Foundation, Grant/Award Number: 1708085QD83; Research Fund for the Doctoral Program of Higher Education of China, Grant/Award Number: 20120141110025; National Natural Science Foundation of China, Grant/Award Numbers: 41074008, 41431069, 41604028; National Science Fund for Distinguished Young Scholars, Grant/Award Number: 41525014

- Received April 1, 2019.
- Revision received December 17, 2019.
- Accepted March 23, 2020.

- © 2020 Institute of Navigation

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.