BDS/GPS Dual Systems Positioning Based on the Modified SR-UKF Algorithm.

The Global Navigation Satellite System can provide all-day three-dimensional position and speed information. Currently, only using the single navigation system cannot satisfy the requirements of the system's reliability and integrity. In order to improve the reliability and stability of the satellite navigation system, the positioning method by BDS and GPS navigation system is presented, the measurement model and the state model are described. Furthermore, the modified square-root Unscented Kalman Filter (SR-UKF) algorithm is employed in BDS and GPS conditions, and analysis of single system/multi-system positioning has been carried out, respectively. The experimental results are compared with the traditional estimation results, which show that the proposed method can perform highly-precise positioning. Especially when the number of satellites is not adequate enough, the proposed method combine BDS and GPS systems to achieve a higher positioning precision.

1. Introduction

With the development of space information technology, some countries are constructing Global Navigation Satellite Systems (GNSS). Now, in addition to USA’s GPS and Russia’s GLONASS, Europe’s Galileo and China’s BeiDou navigation satellite System (BDS) are being built. Japan, India and other countries are also planning to build their own regional navigation satellite systems. Even though it is GPS that is the most developed navigation satellite system and it has many advantages, it also has some disadvantages of system reliability that cannot satisfy the requirements of a single navigation system in certain situations [1]. Recently, the idea of multi-navigation positioning that consists of GPS, GLONASS, Galileo and regional satellite positioning system is gradually getting more interest in the field of satellite navigation. Especially, the combination of GPS and BDS can overcome the deficiency of the single system, and shows better effects on system performance [2].

The most conventional positioning estimation method is the iterative least square method (ILS). Furthermore, extended Kalman filter (EKF) and unscented Kalman filter (UKF) are also used to estimate the positioning data. ILS can solve three-dimensional positioning only when it receives the signals from at least four satellites. This method is simple, and its computing speed is fast, but it has a large linearization error and a low positioning estimation precision. The EKF can only be accurate for a first order Taylor series. There may be a larger nonlinear error, and it needs to compute the Jacobian matrix, in addition to the calculation being difficult and one of the main sources of error [3]. The UKF represents statistical properties of the system by deterministic sampling and avoids the disadvantage that the EKF must compute the Jacobian matrix. Theory shows that the EKF predicts the means correctly up to the second order of Taylor series and covariances up to fourth order. In contrast, the UKF predicts the means and covariances correctly up to the fourth order [4,5].

Currently, unscented Kalman filter (UKF) and square root UKF (SR-UKF) are the widely used nonlinear filtering strategies, their applications are proposed. In the process of filtering, the calculation error exists. The accumulation of the calculation error reduces the filtering precision, and even it can make the error covariance matrices gradually lose their positive semidefiniteness [6]. In order to improve the numerical performance, the SR-UKF was proposed [7]. Here, Cholesky factors of the covariance matrices are directly used to calculate the state estimate, the QR decomposition, as well as Cholesky factor updates. Thus, by this way, the numerical stability can be improved and also the positive semi-definiteness of covariance matrices can usually be guaranteed [8].

In practical applications, the statistic models of the extended noises are difficult to build. In this case, Kalman filter will be invalid, so the adaptive filtering and robust filtering theory is brought and developed [9].

Because the norm of estimation error covariance decreases progressively in the process of filtering, the effects of new observations data for improving state estimation will be weakened. In fact, the changes of system dynamic model are difficult to fully know in advance. As the recent observation data contains more information about the changed system model, a modified SR-UKF algorithm is proposed in order to increase the weight of the new measurement data.

This paper proposes a BDS-GPS system model, using the modified SR-UKF algorithm to perform position estimation, and the experimental results illustrate the effectiveness of our proposed method.

Section 2 lists the general concept of the GNSS positioning. Section 3 introduces the modified SR-UKF algorithm which will be used in our proposed models for GNSS position estimation. Section 4 describes the unification of the reference systems and the time systems of GPS and BDS. Section 5 addresses the developed nonlinear model and the filter implementation . Section 6 includes recent experimental results and provides the comparison of these results from GPS and BDS with ILS, UKF, SR-UKF and the proposed method. Finally, Section 7 summarizes this paper, and puts forward the future developments.

2. GNSS Positioning Overview

GNSS is a worldwide all-weather navigation system which can provide tridimensional position, velocity, and time synchronization to the UTC scale. GNSS considers the earth’s center as the reference point, to determine the position of the receiver antenna in the reference coordinate system. Since the positioning operation requires only one receiver, it is called standalone positioning. The basic principle of the standalone GNSS positioning is taking the observed distance between the GNSS satellite and the user receiver antenna as the benchmark, which is based on the known instantaneous satellites’ coordinates, to determine the position of the corresponding user receiver antenna. According to the different positions of the user receiver antennas, GNSS positioning can be divided into dynamic positioning and static positioning. Currently, GNSS positioning has a variety of modes, such as precise point positioning (PPP) and relative positioning. Precise point positioning uses the precise ephemeris and satellite clock bias data provided from the International GNSS Service (IGS) and calculates the precise user coordinates from the corrected carrier phase and pseudorange. The relative positioning uses single or double difference of the carrier phase, and calculates the relative coordinates comparing to one or several base stations. These technologies can perform high precise positioning, but they need some accessories besides the user GNSS receiver, such as radio, network equipment, and the other GNSS receivers. In a word, they are not pure standalone positioning, due to their high cost and complex setup, and they are unsuited to the daily applications such as car navigation.

GNSS positioning is based on the one-way ranging technique: the propagation time to transmit from satellite to user receiver is measured and multiplied by the signal propagation velocity to obtain satellite-to-user range. The offset of the receiver clock relative to the system time scale should be estimated to position. The measured range between receiver and satellite is referred to as pseudorange, and can be represented as follows: where is the ith satellite’s pseudorange measurement, is the geometric distance receiver-satellite, is the receiver clock offset (scaled by speed of light c ), and contains the residual errors after satellite-based and atmospheric error corrections [10].

Equation (1) is applicable to the single GNSS (i.e., BDS or GPS only), it contains the time scale of the considered system. However, for the multiple constellation case, another unknown variable used to represent the inter-system offset should be further estimated.

3. The Modified SR-UKF Algorithm

The UKF algorithm is the minimum variance estimation based on UT (Unscented Transform). It was first proposed by Julier et al. [11] in 1995. The state distribution here is represented by a number of appropriately chosen points, which is different from the Gaussian Random Variables (GRV) in UKF with deterministic samplings. These points evolve according to the dynamics of the true nonlinear system. Hence, compared with the EKF, the UKF not only has the possibility to improve the estimate precision, but also is easier to be implemented. Moreover, different from EKF, the evaluations of the Jacobian and any order of partial derivatives are not needed in the UKF. Some papers proposed the new EKF and UKF algorithm [12,13,14,15], and were used in GPS positioning [16,17,18,19,20]. The fuzzy adaptive UKF algorithm was applied in spacecraft celestial navigation [21]. When no more than four satellites can be received, a precision of data processing can be obtained by considering the UKF algorithm’s small linearization error.

The UKF is mainly used for an arbitrary nonlinear system, and numerical instability often causes the covariance matrix P to lose its positive definiteness during the filtering procedure. Consequently, the sigma points cannot be correctly calculated, where is a priori estimate of state, is a priori covariance matrix of state, (for L and λ, see Equation (5)). Moreover, in the UKF design, the demanding operation is the evaluation of the square root of the covariance matrices at each time instant for the updated set of sigma points. To solve this problem, the SR-UKF was proposed [7].

Meanwhile, in the application of positioning, the exact knowledge of the noise matrix which is required in the framework of the Kalman filter is usually unknown and time-varying in practice. The inappropriate prior statistics in the Kalman filter cause large estimation errors or even errors possibly diverging. Because of the uncertain process noise, the standard UKF yields poor performance in robustness and tracking accuracy.

In the process of standard filtering, the norm of the estimation error covariance matrix is reduced with time, thus the effects of observations for correction of the state estimation are more and more weakened. As the recent observations contain more information on the changed dynamic system model, in the process, the effects of new observations for the state estimation error must be enhanced, and the effects of the old observations must be reduced.

First, assume that state and measurement equations of the system are discrete time nonlinear systems: where is state vector, is measurement vector, and is zero-mean independent Gaussian white noise, of which the covariance matrix is Q. is zero-mean independent Gaussian white noise of the measurement, of which the covariance matrix is R.

The UKF and SR-UKF algorithms are given in Table 1 and Table 2.

Table 1

Unscented Kalman filter.

Calculate sigma points	χ_t-1=[x^_t-1x^_t-1±(L+λ)P_t-1]
Time-update equations	χt|t-1*=f(χ_t-1)
	x^_t|t-1=∑i=02LWi(m)χi,t|t-1*
	P_t|t-1=∑i=02LWi(c)[χi,t|t-1*-x^t|t-1-]·[χi,t|t-1*-x^t|t-1-]^T+Q
	z_t|t-1=h(χt|t-1*)
	z^_t|t-1=∑i=02LWi(m)z_i,t|t-1
Measurement-update equations	P_z_tz_t=∑i=02LWi(c)[z_i,t|t-1-z^_t|t-1]·[z_i,t|t-1-z^_t|t-1]^T+R
	P_x_tz_t=∑i=02LWi(c)[χi,t|t-1*-x^_t|t-1]·[z_i,t|t-1-z^_t|t-1]^T
	K_t=P_x_tz_tPz_tz_t-1
	x^_t=x^_t|t-1+K_t(z_t-z^_t|t-1)
	P_t=P_t|t-1-K_tP_z_tz_tKtT

Table 2

Square-root unscented Kalman filter.

Calculate sigma points	χ_t-1=[x^_t-1x^_t-1±(L+λ)S_t-1]
Time-update equations	χt|t-1*=f(χ_t-1)
	x^_t|t-1=∑i=02LWi(m)χi,t|t-1*
	S_t|t-1=qr{W1(c)(χ1:2n,t|t-1*-x^_t|t-1)Q}
	S_t|t-1=cholupdate{S_t|t-1,χ0,t|t-1*-x^_t|t-1,W0(c)}
	z_t|t-1=h(χt|t-1*)
	z^_t|t-1=∑i=02LWi(m)z_i,t|t-1
Measurement-update equations	S_z_t=qr{W1(c)(z_1:2n,t|t-1-z^_t|t-1)R}
	S_z_t=cholupdate{S_z_t,z_0,t|t-1-z^_t|t-1,W0c}
	P_x_tz_t=∑i=02LWi(c)[χi,t|t-1*-x^_t|t-1]·[z_i,t|t-1-z^_t|t-1]^T
	K_t=(P_x_tz_t/Sz_tT)/S_z_t
	x^_t=x^_t|t-1+K_t(z_t-z^_t|t-1)
	S_t=cholupdate{S_t|t-1,K_tS_z_t,-1}

Details about the modified SR-UKF algorithm are described as follows:

Initialize with where denotes Cholesky factorization. For a positive definite matrix A, if a matrix X is lower triangularly denoted by , then X is the Cholesky factor of A. The shorthand notation denotes a Cholesky factorization, namely, .

Calculate the sigma points:

Calculate weight coefficient: where is a scaling parameter. The constant ξ determines the spread of the sigma points around a mean of state x, and is usually set to a small positive value (e.g., ). The constant k is a secondary scaling parameter, which is usually set to , and η is used to incorporate prior knowledge of the distribution of x (for Gaussian distributions, is optimal) [5,11].

Time-update equations: where denotes QR decomposition. A matrix is given by , and is orthogonal, and is upper triangular. is the upper triangular part of R. We use the shorthand notation to denote QR decomposition, namely, , is the update to Cholesky factorization, returns the Cholesky factor of , where R is the original Cholesky factorization of A.

Measurement-update equations: where, and S is selected by experience. If S is too large, it will cause filter oscillation. Thus, generally, S is chosen to be slightly greater than one. (In this experiment, .)

4. The Unification of Time and Reference Systems of BDS and GPS

4.1. The Unification of the Two Reference Systems

For the reference systems of BDS’s CGCS2000 and GPS’s WGS84, the definitions are essentially the same [22], and the reference ellipsoids are very similar. The extremely small difference is mainly on the flat rate. Due to this difference, the maximum deviations in the latitude and the altitude are 0.105 mm, respectively, and the maximum deviation of gravity is 0.016 × 10 m/s. These deviations can be neglected under the current accuracy of measurement level. Thus, the GPS and BDS’s position data can be used directly without coordinate transformation.

4.2. The Unification of the Two Time Systems

There is no leap second problem since both BDS and GPS use atomic time. The two systems differ by 1356 weeks and 14 s due to the difference on the starting point. It gives the synchronization parameters between BDS and GPS in the navigation messages. However, the implementation content is unpublished. Thus, the two time systems cannot be fully synchronized directly. Alternatively, we choose a crude way by simply adding 14 s to BDT’s time, and adding an unknown variable to represent the time deviation in different systems [23].

5. The Positioning of BDS and GPS Using the Modified SR-UKF

Nonlinear Kalman filter can be well applied in the GNSS positioning estimation because of its characteristics in which the current state parameter is updated according to the observed value using the predictive value. The system model consists of the process model and the measurement model.

5.1. Process Model

The state model includes the receiver position and velocity coordinates in CGCS2000 coordinate system, and the receiver clock bias, which is related with states and the clock drift caused by the Doppler deviation. It also includes the non-white error in each satellite channel. Thus, the overall system state has 10 fundamental states plus one shaping state for each observable channel: where is the receiver’s position coordinate in CGCS2000, is the receiver velocity coordinate, are clock bias and clock drift bias, respectively, is the clock deviation between GPS and BDS, is the offset between the Doppler shift and pseudorange’s rate, is the non-white error in each satellite channel. Because the error in each satellite is independent, it can be modeled as a first-order Gaussian–Markov process.

Since there is a deviation in the two systems’ times, in order to eliminate positioning error caused by time deviation between BDS and GPS, a variable should be added.

Considering the above state model, and a generic kinematics model for the receiver coordinates, we obtain the associated system model:

The state transition matrix F is given by where A is a time invariant matrix and B is a diagonal matrix, which are given respectively by: where .

is the system driving noise, which is given by where and are noises due to the receiver acceleration and other system interferences, is the driving noise of the ith shaping filters of the channel, and are driving noise for the clock bias model, is the noise of the clock deviation between BDS and GPS, and is the noise of . For each channel i of the received satellite, the non-white component could be modeled using the first-order Gaussian–Markov process since system errors are independent.

The noise matrix C is given by where D is a time invariant matrix and E is a diagonal matrix, which are given respectively by where .

The correspondent process noise covariance matrix is given by with , where is the variances associated with the , is the variances associated with the , and are the process noise variances associated with and , are the variances associated with the shaping filters driving noises and, is defined by following equations: with where [3,17].

5.2. Measurement Model

The pseudoranges and Doppler shifts form the measurements set, and the measurement equations of the BDS’s pseudorange and GPS’s pseudorange are as follows: where are the receiver position coordinates, is the ith BDS satellite’s coordinates, is the jth GPS satellite’s coordinates, is the non-white error in satellite channel i, and is the measurement noise of channel i.

Doppler shifts give information related to the receiver velocity. Doppler is also used in our formulation, modeling it as: where are velocity coordinates of receiver at time t; and is a velocity coordinate of satellite i at time t .

Thus, the system measurement for n satellites is given by: where , , and are the number of measured BDS and GPS satellites.

The estimation flow using the nonlinear filter is shown in Figure 1.

Figure 1

Flow of the position estimation using the nonlinear filter.

6. Experiment and Analysis

In the experiment, a set of BDS/GPS data collected by the OEM-615 (NovAtel Inc., Calgary, AB, Canada) receiver are used. This experiment was carried on School of Electronic Information and Electric Engineering of Shanghai JiaoTong University. The experimental data were collected between 7:30 a.m.–8:10 a.m., 12 October 2015. The sampling period is 1 s.

For the analysis, static receiver data were collected. The mean value of receiver’s 24 h static position data was set as the benchmark. The coordinate of the benchmark in ECEF is .

The number of visible satellites during the experiment’s period is shown in Figure 2. Since the observation environment was good in the outdoor with a clear view of the sky, the number of BDS satellites was stable at 10, GPS satellite numbers changed between six and seven, especially the changes in 600–800 s and 1400–1600 s were more variable. The position estimation is performed using MATLAB in the PC. The raw data was processed using four algorithms: ILS, UKF, SR-UKF and the modified SR-UKF.

Figure 2

The number of visible satellites.

The position estimation errors are shown in Figure 3, Figure 4 and Figure 5.

Figure 3

Position estimation error using BDS’s data.

Figure 4

Position estimation error using GPS’s data.

Figure 5

Position estimation error using BDS/GPS’s data.

The current BDS satellites are distributed in geostationary orbit (GEO) and inclined geosynchronous satellite orbit (IGSO) over China. That is why the BDS satellite number that is received is stable at the present stage. Due to the change in satellite numbers, the positioning results of only using GPS system experiments are not better than BDS single system and BDS/GPS dual systems.

The correspondent root mean square errors (RMSE) are compared in Table 3, Table 4 and Table 5. It is not difficult to find whether BDS, GPS or BDS/GPS positioning, the result of the modified UKF has higher accuracy than ILS, UKF and SR-UKF.

Table 3

Position estimation error using BDS’s data.

	RMSE /m
	X	Y	Z	3D
ILS	0.608	0.595	2.159	2.321
UKF	0.435	0.819	1.929	2.141
SR-UKF	0.434	0.721	1.907	2.085
Modified SR-UKF	0.336	0.357	1.894	1.956

Table 4

Position estimation error using GPS’s data.

	RMSE /m
	X	Y	Z	3D
ILS	3.136	2.371	3.666	5.376
UKF	3.029	2.342	3.592	5.251
SR-UKF	3.015	2.301	3.314	5.036
Modified SR-UKF	2.922	2.118	3.267	4.869

Table 5

Position estimation error using BDS/GPS’s data.

	RMSE /m
	X	Y	Z	3D
ILS	2.049	1.384	1.791	3.053
UKF	2.107	1.355	1.551	2.946
SR-UKF	2.106	1.288	1.543	2.911
Modified SR-UKF	2.062	1.268	1.539	2.869

When the number of the visible satellites is large, this method cannot fully represent its superiority. We have eliminated some of the original data so as to simulate the condition of less visible satellites, the case when the numbers of BDS and GPS were in the 2–3 range.

The positioning with BDS and GPS in bad environments is shown as Figure 6.

Figure 6

Positioning with GPS and BDS under bad conditions.

In a single system, where there are less than four visible satellites, the positioning cannot be obtained by the conventional method, but a multi-system positioning model is a good way to deal with this situation. The number of visible satellites after part of them have been removed is shown in Figure 7. The number of the satellites of BDS remains at three, but that of GPS changes in 2–3.

Figure 7

The number of visible satellites after removes.

When the number of satellites is less than four, the single system cannot be calculated; the advantage of multi-system positioning can now be represented. As long as the sum of BDS and GPS satellites is not less than five, the multi-system can be calculated effectively.

In Figure 8, we can see that the number of satellites is going down, ILS algorithm’s positioning precision is getting worse, and the rapid change in the number of satellites causes unsatisfactory results. Kalman filtering overcomes these difficulties well. Reduction in the number of satellites does not have a significant impact on the positioning result.

Figure 8

Estimation error using BDS/GPS’s Data after removal.

The correspondent root mean square errors of the various methods are compared in Table 6. It is not difficult to find that the result of the modified UKF has the highest accuracy. When the number of BDS or GPS satellites is too few to calculate in each single positioning system, we can adopt the multi-system positioning method. We can get the results by using the ILS algorithm, but the error is too large and it does not meet the requirements of practical application. The Kalman filtering algorithm of the multi-system can still maintain high precision, and it has a huge advantage.

Table 6

Positioning estimation error using BDS/GPS’s data after removal.

	RMSE /m
	X	Y	Z	3D
ILS	5.425	15.954	56.836	59.282
UKF	2.164	4.262	7.425	8.829
SR-UKF	1.475	3.642	6.282	7.409
Modified SR-UKF	0.801	1.567	5.731	5.995

7. Conclusions

This paper has presented the results of a study about the modified SR-UKF algorithm for multi system positioning, which was verified by the real BDS/GPS data. The proposed method is suited to the standalone GNSS positioning with low-cost and wide application. The main contributions are summarized as follows:

The new nonlinear positioning models for the two navigation systems were improved. Considering the differences in dual systems, a new set up for the state variables was proposed. Because there is a deviation in the two systems’ time, in order to eliminate positioning error caused by time deviation between BDS and GPS, a variable must be added. To remedy the shortage that various errors are simply taken as the Gaussian white noise in the iterative least square method, a multi-GNSS positioning model based on the nonlinear filtering was designed by establishing a suitable error model in the positioning system model and taking the nonlinear pseudorange and Doppler observation equation as the measurement equation.

After analyzing the characteristics of BDS and GPS, the unity of the coordinate and time systems of the two systems was considered. The BDT was selected as the time standard, and the CGCS2000 as the coordinate standard.

To reduce the bad effect of the old measurement data on the filtering and increase the weight of the new measurement data, the modified SR-UKF algorithm was used. The proposed algorithm gradually decreases the weighting on the old measurement data and increases it on the new measurement data, correspondingly, which overcomes the filter divergence effectively.

The new model and method were processed by ILS, UKF, SR-UKF and the proposed method. The experimental results show that the BDS/GPS systems positioning estimated by the proposed algorithm performs the best.

The proposed method also can be used in multi-system positioning in urban canyon environments. In a single system when there are less than four visible satellites, the positioning cannot be obtained by the conventional method. However, using the proposed nonlinear positioning models, the high precision positioning can be solved as long as the sum of the satellites is no less than five in the dual system. In addition, the positioning accuracy is much higher than the result obtained by the iterative least square algorithm.

Future works include applying the new model and new method to a dynamic environment, adding more GNSS systems into the environment, and making it more stable under the bad positioning conditions in which the number of satellites is not enough. In order to achieve the balance between positioning accuracy and computational complexity, we will investigate the satellite selection problem.