A Window-Recursive Approach for GNSS Kinematic Navigation Using Pseudorange and Doppler Measurements

Abstract

In kinematic Global Navigation Satellite Systems (GNSS) navigation, the Kalman Filter (KF) solution relies, to a great extent, on the quality of the dynamic model that describes the moving object's motion behaviour. However, it is rather difficult to establish a precise dynamic model that only connects the previous state and the current state, since these high-order quantities are usually unavailable in GNSS navigation receivers. To overcome such limitations, the Window-Recursive Approach (WRA) that employs the previous multiple states to predict the current one was developed in Zhou et al., (2010). Its essence is to adaptively fit the moving object's motion behaviour using the multiple historical states in a short time span. Up to now, the WRA method has been performed only using GNSS pseudorange measurements. However, in GNSS navigation fields, the strength of pseudorange observation model is usually weak due to various reasons, e.g., multi-path delay, outliers, insufficient visible satellites. As an important complementary measurement, Doppler can be used to aid Position and Velocity (PV) estimation. In this contribution, implementation of WRA will be developed using the pseudorange and Doppler measurements. Its corresponding state transition matrix is constructed based on the Newton's Forward Difference Extrapolation (NFDE) and Definite Integral (DI) methods for the efficient computation. The new implementation of WRA is evaluated using the real kinematic vehicular GNSS data with two sampling rates. The results show that:

1. (i) aided by GNSS Doppler measurement, the new implementation of WRA significantly improves the accuracy compared with the pseudorange-only WRA.

2. (ii) In high sampling rate, the WRA works best in the case of 2 epochs in time window, while in the low sampling rate, it obtains better solutions if more epochs involved in time window.

3. (iii) Compared with KF with constant velocity dynamic model, the WRA demonstrates better in the self-adaptation and validity.

4. (iv) As a benefit of WRA itself, the NFDE/DI-based state transition matrix for WRA can be previously computed offline without increasing the computation burdens.

1. INTRODUCTION

The Global Navigation Satellite Systems (GNSS) can provide real-time Position and Velocity (PV) information suitable for most land, marine, and aircraft navigation applications (Parkinson and Spilker 1996). In GNSS kinematic positioning and navigation, the pseudorange measurement is preferred to the carrier-phase measurement due to its exemption of unknown integer cycle. In general, the metre and sub-metre level kinematic solutions are achievable from the single point positioning and the differential GNSS technique, in real time respectively. As an effective and important complement to GNSS kinematic navigation, the GNSS Doppler measurements provide the velocity information of moving objects (Parkinson and Spilker 1996; Bruton et al., 1999; Zhang et al., 2006; Remondi, 2004; Zhou et al., 2011). In GNSS kinematic navigation and positioning, it still faces many challenges due to, for instance, the signal degradation, the blockages and the colored noise etc. Extensive research efforts have been done towards optimizing the solution and maximizing the reliability in past two decades (Teunissen 1990; Wang 2000; Hewitson and Wang 2007; Kuusniemi et al., 2004).

As an efficient implementation of the sequential least squares adjustment with relative constraints between-state, the Kalman Filter (KF) is extensively applied in the navigation computations assuming that the observational error and the prediction error are uncorrelated and both normally distributed with zero means (Brown and Hwang 1997; Simon, 2006). In real applications, the unexpected errors inevitably exist in either the dynamic model or the observation model. Therefore, to balance the contributions of the observation and the prediction to the solution, the adaptive and robust filtering techniques were developed based on the robust M-estimation, optimal theory and Bayesian estimation (Koch and Yang 1998; Yang et al., 2001, 2003, 2005; Wang 2000, Hewitson and Wang 2007; Geng and Wang 2008).

The KF approach relies, to a great extent, on the quality of dynamic model that describes the moving object's motion behaviour. In the existing literatures, the Constant Velocity (CV) and the Constant Acceleration (CA) dynamic models are popularly employed in navigation. However, the CV model is too simple to describe the motion detail and then suits only for those low-dynamic or high sampling rate systems. Although the CA model can capture the motion in more details, the acceleration cannot be directly obtained from the GNSS receivers. Furthermore, both two models would be problematic to describe the moving object's manoeuvring that is a complicated and unpredicted dynamic process. Thorp (1973) proposed to model this dynamic process of manoeuvring as a stochastic departure from a nominated base trajectory to obtain the reasonable performances, while Singer (1970) solved for such dynamic process as a noisy component by using a known exponential auto-correlation acceleration function. In addition, an extended ‘Singer’ model was then developed permitting the time-correlation in both the velocity and the acceleration components (Moore and Wang 2003). Another dynamic model based on ‘kinematic’ method was applied in KF, which introduces the variation of moving object's positions to reflect the mean velocity in two consecutive epochs (Schwarz et al., 1989).

These dynamic models aforementioned are established based on state of the previous one epoch only, which is impossible to provide the sufficient information for describing the moving object's motion characteristic of current epoch, thus leading to the increased degradation and instability risks in kinematic solutions. Therefore, it is necessary to include more historical information in the dynamic model. The Windowing-Recursive Approach (WRA) provides a new strategy to estimate the state of vehicles (Zhou et al., 2010). Its basic principle is that the trajectory of a vehicle can be reasonably fitted by using the positions of latest several epochs in a short time span. Therefore the position at the current epoch can be predicted more reliably with the positions of the multiple historical epochs rather than the latest one epoch. In WRA, the dynamic model of vehicle is adaptively recovered in real time without any assumption on the dynamic characteristics, like CV model. The WRA is flexible and can reduce to KF if the window consists only of one epoch. However, in the first WRA literature of Zhou et al., (2010), the WRA was implemented with the GNSS pseudorange measurements and the efficient state transition matrix was constructed only for the case where the state vector consists only of three-dimensional position. Nevertheless, velocity needs to be estimated in real time besides position in most of GNSS navigation applications, such as aircrafts and missiles. Therefore, the implementation of WRA should be extended from the three-dimensional position state to the six-dimensional PV state to estimate the reliable and accurate real-time solutions by taking full advantages of velocity information. In this contribution, the implementation of WRA with six-dimensional PV state vector will be developed using the pseudorange and Doppler measurements. Its corresponding state transition matrix is constructed based on the Newton's Forward Difference Extrapolation (NFDE) and Definite Integral (DI) methods for the efficient computation.

The rest of paper is organized as follows. In Section 2, Kalman filtering and the WRA are introduced in brief. In Section 3, the implementation of WRA is developed based on the pseudoranges and Doppler measurements. In Section 4, we will construct the adaptive, efficient and practical state transition matrices for the proposed WRA implementation based on the NFDE and DI numerical methods. A real vehicular GNSS experiment is introduced in Section 5 to demonstrate the performance of the new WRA implementation. Finally, some concluding remarks are given in Section 6.

2. KALMAN FILTERING AND WINDOW-RECURSIVE APPROACH

2.1. Kalman Filtering

As an efficient realization of the sequential least-squares adjustment, the KF has been widely used in the GNSS navigation computations. A linear dynamic model and an observation model are involved in the KF:

(1a)$${\bi x}_k = {\bi \Phi} _{k,k - 1} {\bi x}_{k - 1} + {\bi w}_k $$

(1b)$${\bi l}_k = {\bi H}_k {\bi x}_k + {\bi \varepsilon} _k $$

where:

The subscripts denote the epoch number.

x is the state vector to be estimated.

l is the measurement vector.

w and ε are the process noise in dynamic model and the measurement noise with the zero means and the covariance matrices of Σ _w and R, respectively.

Φ is the state transition matrix that transfers the state vector from the (k − 1)th epoch to the kth epoch.

H is the design matrix connecting the state vector with the observation vector.

The sequential formulae in KF are given:

(2a)$${\bi \bar x}_k = {\bi \Phi} _{k,k - 1} {\bi \hat x}_{k - 1} $$

(2b)$${\bi \Sigma} _{{\bi \bar x}_k} = {\bi \Phi} _{k,k - 1} {\bi \Sigma} _{{\bi \hat x}_{k - 1}} {\bi \Phi} _{K,k - 1}^T + {\bi \Sigma}_{{\bi w}_{k}} $$

(2c)$${\bi K}_k = {\bi \Sigma} _{{\bi \bar x}_k} {\bi H}_k^T \left( {{\bi H}_k {\bi \Sigma} _{{\bi \bar x}_k} {\bi H}_k^T + {\bi R}_k} \right)^{ - 1} $$

(2d)$${\bi \hat x}_k = {\bi \bar x}_k + {\bi K}_k \left( {{\bi l}_k - {\bi H}_k {\bi \bar x}_k} \right)$$

(2e)$${\bi \Sigma} _{{\bi \hat x}_k} = \left( {{\bi I} - {\bi K}_k {\bi H}_k} \right){\bi \Sigma} _{{\bi \bar x}_k} $$

where:

I is the identity matrix with same dimensions of ${\bi \Sigma} _{{\bi \bar x}_k} $.

${\bi \bar x}_k $${\bi \Sigma} _{{\bi \bar x}_k} $ denotes the predicted state vector and its covariance matrix, respectively.

K_k is the gain matrix. ${\bi \hat x}_k $ and ${\bi \Sigma} _{{\bi \hat x}_k} $ are the posterior KF estimate and its covariance matrix.

If the dynamic model or/and observation model is nonlinear, the Extended Kalman Filtering (EKF) would be applied where these nonlinear models are linearized (Brown and Hwang 1997; Simon 2006). The EKF can be viewed as the conventional KF in the first-order approximated linear model. These approximations, however, may introduce the large errors in the true posterior mean and covariance of the transformed Gaussian Random Variable (GRV), resulting in the sub-optimal performance or even filtering divergence. To overcome the pitfalls of EKF, the sigma-points and Unscented Kalman Filtering (UKF) approaches were developed (Julier and Uhlmann 1997; Julier et al., 2000; Jwo and Lai 2009) and have been successfully applied in target tracking field. The state distribution is again represented by a GRV but specified by using a minimal set of carefully chosen sample points. These sampled points can adequately capture posterior mean and covariance accurately to the 3^rd order (Taylor series expansion) for any nonlinearity when propagated through the true non-linear system. Since the nonlinearity of the GNSS observational model is so weak that we can basically neglect its influence on estimation accuracy, the EKF is preferable in GNSS navigation (Jwo and Lai 2008).

2.2 A Brief Review on Window-Recursive Approach

In traditional navigation applications, the dynamic model is established based on the state of the previous one epoch only, which is impossible to provide the sufficient information for describing the moving object's motion characteristics of current epoch, thus leading to the increased degradation and instability risks in kinematic solutions. To overcome this limitation, the WRA was proposed in Zhou et al., (2010). Its principle is that the position at the current epoch is predicted more reliably with the positions of the multiple historical epochs than with the position only of the latest epoch. Assuming that the time-window length contains n epochs, the dynamic model is constructed as:

(3)$${\bi x}_k = {\bi \Phi} _{\left( {k,k - n:k - 1} \right)} {\bi x}_{\left( {k - n:k - 1} \right)} + {\bi w}_k $$

where

x_{(k − n:k − 1)} is the vector consisting of the stacked state vectors from the epoch (k − n) to the epoch (k − 1).

Φ_{(k ,k − n:k − 1)} is the transition matrix that transfers the state information of the previous n epochs into that of the current one.

Let ${\bi \Sigma} _{{\bi \hat x}_{\left( {k - n:k - 1} \right)}} $ be the covariance matrix of the historical state vector ${\bi \hat x}_{\left( {k - n:k - 1} \right)} $, the predicted position vector at epoch k and its corresponding covariance are derived:

(4)$${\bi \bar x}_k = {\bi \Phi} _{\left( {k,k - n:k - 1} \right)} {\bi \hat x}_{\left( {k - n:k - 1} \right)} $$

and:

(5)$${\bi \Sigma} _{{\bi \bar x}_k} = {\bi \Phi} _{\left( {k,k - n:k - 1} \right)} {\bi \Sigma} _{{\bi \hat x}_{\left( {k - n:k - 1} \right)}} {\bi \Phi} _{\left( {k,k - n:k - 1} \right)}^T + {\bi \Sigma} _{{\bi w}_k} $$

respectively.

In the linear(ised) observation Equation (1b) at epoch k, the covariance matrix of measurements is expressed as R_k=σ ₀²P_k⁻¹. Here, σ ₀² is a known prior variance scalar and P_k is a weight matrix. By minimizing the weighted quadratic measurement residual vector ${\bi v}_k = {\bi H}_k {\bi \hat x}_k - {\bi l}_k $ and the weighted quadratic predicted residual vector ${\bi \bar v}_{{\bi x}_k} = {\bi \hat x}_k - {\bi \bar x}_k $, the Bayesian risk function is established by:

(6)$$\mathop {\min} \limits_{{\bi x}_k} :{\bi v}_k^T {\bi P}_k {\bi v}_k + {\bi \bar v}_{{\bi \bar x}_k} ^T {\bi P}_{{\bi \bar x}_k} {\bi \bar v}_{{\bi \bar x}_k} $$

where ${\bi P}_{{\bi \bar x}_k} = \sigma _0^2 {\bi \Sigma} _{{\bi \bar x}_k} ^{ - 1} $ denotes the weight matrix of the predicted state vector.

The solution of minimization problem Equation (6) is:

(7)$${\bi \hat x}_k = {\bi \bar x}_k + {\bi K}_k \left( {{\bi l}_k - {\bi H}_k {\bi \bar x}_k} \right)$$

(8)$${\bi \Sigma} _{{\bi \hat x}_k} = \left( {{\bi I} - {\bi K}_k {\bi H}_k} \right){\bi \Sigma} _{{\bi \bar x}_k} $$

(9)$${\bi K}_k = {\bi \Sigma} _{{\bi \bar x}_k} {\bi H}_k^T \left( {{\bi H}_k {\bi \Sigma} _{{\bi \bar x}_k} {\bi H}_k^T + {\bi R}_k} \right)^{ - 1} $$

The correlation between ${\bi \hat x}_k $ and ${\bi \hat x}_{\left( {k - n:k - 1} \right)} $ should be rigorously processed when the time window moves forward. Hereby, the covariance matrix ${\bi \Sigma} _{{\bi \hat x}_{\left( {k - n:k} \right)}} $ is expressed as:

(10)$${\bi \Sigma} _{{\bi \hat x}_{\left( {k - n:k} \right)}} = \left[ {\matrix{ {{\bi \Sigma} _{{\bi \hat x}_{\left( {k - n:k - 1} \right)}}} & {{\bi \Sigma} _{_{{\bi \hat x}_{_k} {\bi \hat x}_{\left( {k - n:k - 1} \right)}}} ^T} \cr {{\bi \Sigma} _{{\bi \hat x}_k {\bi \hat x}_{\left( {k - n:k - 1} \right)}}} & {{\bi \Sigma} _{{\bi \hat x}_k}} \cr}} \right]$$

Inserting Equation (4) into Equation (7) yields:

(11)$${\bi \hat x}_k = \left( {{\bi I} - {\bi K}_k {\bi H}_k} \right){\bi \Phi} _{\left( {k,k - n:k - 1} \right)} {\bi \hat x}_{\left( {k - n:k - 1} \right)} + {\bi K}_k {\bi l}_k $$

The covariance matrix of ${\bi \hat x}_k $ and ${\bi \hat x}_{\left( {k - n:k - 1} \right)} $ is derived from Equation (12) in terms of the error propagation law as:

(12)$${\bi \Sigma} _{{\bi \hat x}_k {\bi \hat x}_{\left( {k - n:k - 1} \right)}} = \left( {{\bi I} - {\bi K}_k {\bi H}_k} \right){\bi \Phi} _{\left( {k,k - n:k - 1} \right)} {\bi \Sigma} _{{\bi \hat x}_{\left( {k - n:k - 1} \right)}} $$

With the window moving, the new epoch will be introduced and the first epoch of window will be removed. It is rather easy to update the window vector ${\bi \hat x}_{\left( {k - n + 1:k} \right)} $ and its corresponding covariance matrix ${\bi \Sigma} _{{\bi \hat x}_{\left( {k - n + 1:k} \right)}} $ which is a sub-matrix of Equation (10). One issue should be pointed out that the WRA can reduce to KF when window length n = 1. The WRA implementation from epoch (k − 1) to k is presented in Figure 1.

Figure 1. The computation flowchart of WRA.

3. IMPLEMENTATION OF WRA BASED ON PSEUDORANGE AND DOPPLER MEASUREMENTS

Generally, the GNSS navigation receivers output PV quantities in real time. By utilising the PV solutions generated from GNSS receivers and the Equation (3), the WRA can be directly implemented based on Equations (4) ∼ (12). However, when the number of tracked satellites is less than four, the GNSS receivers’ built-in software cannot generate epoch-wise PV solutions anymore. Nevertheless, the raw GNSS measurements (e.g., pseudorange and Doppler) can still contribute to the PV estimation even if the number of satellites is less than four (Han and Wang 2012). Moreover, the raw measurements are more efficient and fundamental in quality control of kinematic system, since raw measurements can be used to identify and correct the errors of either kinematic model or observational model based on robust theory in quality control of navigation system (Yang et al., 2004).

The GNSS pseudorange observable equation is symbolized at epoch t:

(13)$$O_r^s \left( t \right) = \left\| {{\bi p}_r \left( t \right) - {\bi p}^s \left( t \right)} \right\| + c\left( {{\rm\Delta} t_r \left( t \right) - {\rm\Delta} t^s \left( t \right)} \right) + {\rm\Delta} Ion_r^s \left( t \right) + {\rm\Delta} Trop_r^s \left( t \right) + {\rm\Delta} Mp_r^s \left( t \right) + \varepsilon _{O_r^s} \left( t \right)$$

where:

The superscript s denotes the satellite and the subscript r the receiver.

O is the pseudorange observable.

‖p_r−p^s‖ is the satellite-to-receiver distance from satellite s to receiver r.

Δt _r and Δt^s are the receiver and satellite clock offsets with respect to the Global Positioning System Time (GPST) in case of GPS system/constellation, respectively.

c is the speed of light in vacuum.

ΔIon and ΔTrop reflect signal propagation delays due to the ionosphere and the troposphere, respectively.

ΔMp is the multi-path effect.

$\varepsilon_{O_r^s} $ is the un-modelled pseudorange observation error.

Differentiating the pseudorange O _r^s with respect to time t yields the Doppler observable:

(14)$$\eqalign{\lambda D_{rk}^s = & \left. {\displaystyle{{dO_r^s \left( t \right)} \over {dt}}} \right|_{t = k} \cr {\rm} = & \displaystyle{d \over {dt}}\left( {\left\| {{\bi p}_r \left( t \right) - {\bi p}^s \left( t \right)} \right\| + c\left( {{\rm\Delta} t_r \left( t \right) - {\rm\Delta} t^s \left( t \right)} \right) + {\rm\Delta} Ion_r^s \left( t \right) + {\rm\Delta} Trop_r^s \left( t \right) + {\rm\Delta} Mp_r^s \left( t \right) + \varepsilon _{O_r^s \left( t \right)}} \right) \right|_{t = k}\cr {\rm} = & \displaystyle{{({\bi p}_{rk} - {\bi p}_k^s})^T \over {\left\| {{\bi p}_{rk} - {\bi p}_k^s} \right\|}}\left( {{\bi v}_{rk} - {\bi v}_k^s} \right) + c\left( {{\rm\Delta} \dot t_{rk} - {\rm\Delta} \dot t_k^s} \right) + {\rm\Delta} I\dot on_{rk}^s + {\rm\Delta} Tr\dot op_{rk}^s + {\rm\Delta} \dot Mp_{rk}^s + \varepsilon _{\dot O_{rk}^s}} $$

where:

D is the Doppler observable.

λ is the wavelength of carrier.

v is the velocity vector.

Pseudorange and Doppler observation Equations (13) and (14) can be written in a matrix form for a single receiver at epoch k as follows:

(15)$${\bi l}_k = h\left( {{\bi x}_k} \right) + {\bi \varepsilon} _k $$

where:

$${\bi l}_k = \left[ {\matrix{ {{\bi O}_k} \cr {\lambda {\bi D}_k} \cr}} \right],\quad {\bi O}_k = \left[ {\matrix{ {O_{rk}^{s_1}} & {O_{rk}^{s_2}} & \ldots & {O_{rk}^{s_{m_O}}} \cr}} \right]^T $$

$${\bi D}_k = \left[ {\matrix{ {D_{rk}^{s_1}} & {D_{rk}^{s_2}} & \ldots & {D_{rk}^{s_{m_{_D}}}} \cr}} \right]^T $$

m _O and m _D are the number of observed pseudorange and Doppler measurements respectively.

ε_k is the noise vector of l.

${\bi x}_k = \left[ {\matrix{ {{\bi p}_{rk}} \cr {{\bi v}_{rk}} \cr}} \right]$, p_rk and v_rk are the position and velocity vectors, respectively.

h(·) represents the nonlinear relationships Equations (13) and (14) between pseudorange and x_k, and between Doppler and x_k, respectively.

Its explicit expression is:

(15a)$$h\left( {{\bi x}_k} \right) = \left[ {\matrix{ {h_O \left( {{\bi x}_k} \right)} \cr {h_D \left( {{\bi x}_k} \right)} \cr}} \right]$$

and:

(15b)$$h_O \left( {{\bi x}_k} \right) = \left[ {\matrix{ {\left\| {{\bi B}_1 {\bi x}_k - {\bi p}_k^{s_1}} \right\| + c\left( {{\rm\Delta} t_{rk} - {\rm\Delta} t_k^{s_1}} \right) + {\rm\Delta} Ion_{rk}^{s_1} + {\rm\Delta} Trop_{rk}^{s_1} + {\rm\Delta} Mp_{rk}^{s_1}} \cr \vdots \cr {\left\| {{\bi B}_1 {\bi x}_k - {\bi p}_k^{s_{m_O}}} \right\| + c\left( {{\rm\Delta} t_{rk} - {\rm\Delta} t_k^{s_{m_O}}} \right) + {\rm\Delta} Ion_{rk}^{s_{m_O}} + {\rm\Delta} Trop_{rk}^{s_{m_O}} + {\rm\Delta} Mp_{rk}^{s_{m_O}}} \cr}} \right]$$

(15c)$$ \eqalign { \tab \hskip -4pt h_D \left( {{\bi x}_k} \right) = \cr \tab \hskip -5pt \left[ {\matrix{ {\displaystyle{{({\bi B}_1 {\bi x}_k - {\bi p}_k^{s_1})^T} \over {\left\| {{\bi B}_1 {\bi x}_k - {\bi p}_k^{s_1}} \right\|}}\left( {{\bi B}_2 {\bi x}_k - {\bi v}_k^{s_1}} \right) + c\left( {{\rm\Delta} \dot t_{rk} - {\rm\Delta} \dot t_k^{s_1}} \right) + {\rm\Delta} I\dot on_{rk}^{s_1} + {\rm\Delta} Tr\dot op_{rk}^{s_1} + {\rm\Delta} \dot Mp_{rk}^{s_1}} \cr \vdots \cr {\displaystyle{{({\bi B}_1 {\bi x}_k - {\bi p}_k^{s_{m_{_D}}})^T} \over {\left\| {{\bi B}_1 {\bi x}_k - {\bi p}_k^{s_{m_{_D}}}} \right\|}}\left( {{\bi B}_2 {\bi x}_k - {\bi v}_k^{s_{m_{_D}}}} \right) + c\left( {{\rm\Delta} \dot t_{rk} - {\rm\Delta} \dot t_k^{s_{m_{_D}}}} \right) + {\rm\Delta} I\dot on_{rk}^{s_{m_{_D}}} + {\rm\Delta} Tr\dot op_{rk}^{s_{m_{_D}}} + {\rm\Delta} \dot Mp_{rk}^{s_{m_{_D}}}} \cr}} \right]$$

with ${\bi B}_1 = \left[ {\matrix{ {{\bi I}_{3 \times 3}} & {{\bi 0}_{3 \times 3}} \cr}} \right]$ and ${\bi B}_2 = \left[ {\matrix{ {{\bi 0}_{3 \times 3}} & {{\bi I}_{3 \times 3}} \cr}} \right]$, respectively.

It should be pointed out that the systematic errors (such as ionospheric and tropospheric delays) in Equations (13) and (14) can be processed in two ways:

1. (i) Correct these systematic errors with their empirical models.

2. (ii) Use the differenced GNSS model to eliminate these errors.

In this paper, (ii) is employed to eliminate the systematic errors. In addition, considering the weak nonlinearity of the GNSS observational model, similar to the EKF, the WRA equations for the GNSS kinematic and observation model can be conducted as follows:

(16)$${\bi \hat x}_k = {\bi \bar x}_k + {\bi K}_k \left[ {{\bi l}_k - h\left( {{\bi \bar x}_k} \right)} \right]$$

(17)$${\bi \Sigma} _{{\bi \hat x}_k} = \left( {{\bi I} - {\bi K}_k {\bi H}_k} \right){\bi \Sigma} _{{\bi \bar x}_k} $$

(18)$${\bi K}_k = {\bi \Sigma} _{{\bi \bar x}_k} {\bi H}_k^T \left( {{\bi H}_k {\bi \Sigma} _{{\bi \bar x}_k} {\bi H}_k^T + {\bi R}_k} \right)^{ - 1} $$

where ${\bi H}_k = \left. {\displaystyle{{\partial h} \over {\partial {\bi x}_k}}} \right|_{\left( {{\bi x}_k = {\bi \bar x}_k} \right)} $ is the Jacobian matrix.

Without affecting the accuracy of Equation (14), $\displaystyle{{{\bi p}_r - {\bi p}^s} \over {\left\| {{\bi p}_r - {\bi p}^s} \right\|}}$ can be replaced by $\displaystyle{{{\bi p}_{rk0} - {\bi p}_k^s} \over {\left\| {{\bi p}_{rk0} - {\bi p}_k^s} \right\|}}$ and p_rk0 is the approximate position vector of receiver, then H_k is readily derived:

(19)$$\eqalign{{{\bi H}_k = \left[ {\matrix{ {{\bi H}_{{\bi O}_k}} & {{\bi 0}_{m_P \times 3}} \cr {{\bi 0}_{m_D \times 3}} & {{\bi H}_{{\bi D}_k}} \cr}} \right],\hskip5pt{\rm} {\bi H}_{{\bi O}_k} = \left[ {\openup -3pt\matrix{ {\left( {\displaystyle{{{\bi p}_{rk0} - {\bi p}_k^{s_1}} \over {\left\| {{\bi p}_{rk0} - {\bi p}_k^s} \right\|}}} \right)^T} \cr {\left( {\displaystyle{{{\bi p}_{rk0} - {\bi p}_k^{s_2}} \over {\left\| {{\bi p}_{rk0} - {\bi p}_k^s} \right\|}}} \right)^T} \cr \ldots \cr {\left( {\displaystyle{{{\bi p}_{rk0} - {\bi p}_k^{s_{m_P}}} \over {\left\| {{\bi p}_{rk0} - {\bi p}_k^{s_{m_P}}} \right\|}}} \right)^T} \cr}} \right],}\cr {{\rm} {\bi H}_{{\bi D}_k} = \left[ {\openup -3pt\matrix{ {\left( {\displaystyle{{{\bi p}_{rk0} - {\bi p}_k^{s_1}} \over {\left\| {{\bi p}_{rk0} - {\bi p}_k^s} \right\|}}} \right)^T} \cr {\left( {\displaystyle{{{\bi p}_{rk0} - {\bi p}_k^{s_2}} \over {\left\| {{\bi p}_{rk0} - {\bi p}_k^s} \right\|}}} \right)^T} \cr \ldots \cr {\left( {\displaystyle{{{\bi p}_{rk0} - {\bi p}_k^{s_{m_D}}} \over {\left\| {{\bi p}_{rk0} - {\bi p}_k^{s_{m_D}}} \right\|}}} \right)^T} \cr}} \right]\hskip3pt\cr}$$

where ${\bi H}_{{\bi O}_k} $ and ${\bi H}_{{\bi D}_k} $ are the coefficient matrices of pseudorange and Doppler measurements respectively.

Because of its first-order linearization, the estimate of WRA is theoretically suboptimal although the second-order term is very small. It is easy to calculate the covariance matrix of ${\bi \hat x}_k $ and ${\bi \hat x}_{\left( {k - n:k - 1} \right)} $ by:

(20)$${\bi \Sigma} _{{\bi \hat x}_k {\bi \hat x}_{\left( {k - n:k - 1} \right)}} = \left( {{\bi I} - {\bi K}_k {\bi H}_k} \right){\bi \Phi} _{\left( {k,k - n:k - 1} \right)} {\bi \Sigma} _{{\bi \hat x}_{\left( {k - n:k - 1} \right)}} $$

4. STATE TRANSITION MATRIX MODELLING

Different from the conventional KF model, the WRA use the information of multiple historical epochs rather than just the latest one in the dynamic model. Thus the state transition matrix is relatively complicated. To improve the computation efficiency, two different state transition matrices have been established based on Newton's Forward Difference Extrapolation (NFDE) and Polynomial Fitting (PF) numerical methods in the case when only three-dimensional position is involved in the state vector (Zhou et al., 2010). Similarly, for the new WRA implementation in this paper, we need of course to construct their corresponding state transition matrices for efficient computations. It is emphasized again that in the new implementation, the state vector is of six-dimensional position and velocity rather than only three-dimensional position.

Here, the establishment of the six-state model will be realized in two steps. First, since the velocity derived from Doppler is commonly more precise than that differentiated by position (e.g., the CV model where the velocity is derived based on the position difference between two consecutive epochs), we employ the velocity estimated by Doppler to predict the current velocity with the NFDE method in a time window. Secondly, we utilize the velocity extrapolation function of time in the window to calculate the position difference between the current epoch and the previous one by Definite Integral (DI) of that velocity function, and then the current position can be reckoned by the position of previous one epoch and the position difference.

4.1. Velocity Derived from the Doppler Based on the NFDE Model

Similar to the extrapolation model in Zhou et al., (2009, 2010), the velocity vector v can be assumed as a function of time t, i.e.,v_i=v(t _i), where v_i denotes the velocity vector at epoch t _i. By neglecting the high-order remainder terms, the velocity at an arbitrary time can be expressed by the n-th order of NFDE model as:

(21)$$\eqalign{{\bi v}\left( {t_{k - n + s}, n} \right) = & {\bi v}\left( {t_{k - n} + s\delta t,n} \right) \cr = & {\rm\Delta} ^1 {\bi F}_0 + \displaystyle{s \over {1!}}{\rm\Delta} ^2 {\bi F}_0 + \displaystyle{{s\left( {s - 1} \right)} \over {2!}}{\rm\Delta} ^3 {\bi F}_0 + \ldots + \displaystyle{{s\left( {s - 1} \right) \ldots \left( {s - n + 2} \right)} \over {n!}}{\rm\Delta} ^n {\bi F}_0} $$

where:

δt denotes the time interval between two consecutive epochs.

s is the s-th extrapolated epoch s ⩾ n in the time window and n is the order of forward extrapolation.

${\rm\Delta} ^m {\bi F}_0 = \sum\limits_{i = 0}^{m - 1} {\left( { - 1} \right)^{m - 1 - i} C_{m - 1}^i {\bi v}_{k - n + i}} $ with m∈1,\2,\…\n, $C_s^n = \displaystyle{{s!} \over {\left( {s - n} \right)!n!}}$ and (·)! is the factorial operator.

Let s=n, then the velocity at epoch k is conducted as (Zhou et al., 2010):

(22)$$\eqalign{{\bi v}_k = & {\bi v}\left( {t_{k - n + s}, n} \right) = {\bi Jv}_{\left( {k - n:k - 1} \right)} = \sum\limits_{\,j = 1}^n {J_{n - j + 1}} {\bi v}_{k - j} \cr J_j = & \sum\limits_{i = j - 1}^{n - 1} {\left( { - 1} \right)^{i - j + 1} C_i^{\,j - 1}} C_n^i \quad j \in \left\{ {1,{\rm} 2, \ldots {\rm,} n} \right\}} $$

where:

${\bi v}_{\left( {k - n:k - 1} \right)} = \left[ {\matrix{ {{\bi v}_{\left( {k - n} \right)}^T} & {{\bi v}_{\left( {k - n + 1} \right)}^T} & \ldots & {{\bi v}_{\left( {k - 1} \right)}^T} \cr}} \right]^T $ represents a 3n column vector of velocities for all n epochs in the window.

${\bi J} = \left[ {\matrix{ {J_1} & {J_2} & \ldots & {J_n} \cr}} \right] \otimes {\bi I}_{3 \times 3} $ is a 3 × 3n constant coefficient matrix with ⊗ being the Kronecker product.

In real applications, the coefficient matrix J is primarily computed offline and its elements for 1≤n ≤ 5 are presented in Table 1.

Table 1. The coefficient matrix J for the different window lengths (n = 1 to 5)

4.2 Position Derived from Velocity Based on the DI Method

As is well-known in mathematics, in general the differential algorithm augments error while the integral reduces error. Thereby, the definite integral of velocity is introduced to diminish the position errors. At epoch k, position vector can be expressed by position of epoch k − 1 and the position difference from epoch k − 1 and k:

(23)$${\bi p}_k = {\bi p}_{k - 1} + {\rm\Delta} {\bi p}_{k,k - 1} $$

where Δp_{k ,k−1} is the position difference vector from epoch k − 1 to k.

It is rather easy to obtain the position difference of two consecutive epochs which can be recovered by the definite integral of velocity from previous one epoch to current epoch:

(24)$${\rm\Delta} {\bi p}_{k,k - 1} = \int_{t_{\left( {k - 1} \right)}} ^{t_k} {{\bi v}\left( {t_{k - n + s}, n} \right)} dt_s = \int_{k - 1}^k {{\bi v}\left( {t_{k - n + s}, n} \right)} ds$$

Inserting Equation (21) into Equation (24), we obtain the linear combination of velocity vectors in the time window:

(25)$${\rm\Delta} {\bi p}_{k,k - 1} = {\bi Gv}_{\left( {k - n:k - 1} \right)} = \sum\limits_{i = 1}^n {G_{n - i + 1}} {\bi v}_{k - i} $$

Similar to J, the coefficient matrix G is also computed offline and its elements for 1≤n ≤ 5 are presented in Table 2.

Table 2. The coefficient matrix G for the different window lengths (n = 1 to 5)

Finally, the following dynamic model is established based on Equations (22) and (23):

(26)$$\left[ {\matrix{ {{\bi p}_k} \cr {{\bi v}_k} \cr}} \right] = \left[ {\matrix{ {{\bi p}_{k - 1} + {\bi Gv}_{\left( {k - n:k - 1} \right)}} \cr {{\bi Jv}_{\left( {k - n:k - 1} \right)}} \cr}} \right] + \left[ {\matrix{ {{\bi w}_{k1}} \cr {{\bi w}_{k2}} \cr}} \right]$$

Let ${\bi x}_k = \left[ {\matrix{ {{\bi p}_k} \cr {{\bi v}_k} \cr}} \right]$, ${\bi x}_{\left( {k - n:k - 1} \right)} = \left[ {\matrix{ {{\bi x}_{k - n}^T} & {{\bi x}_{k - n + 1}^T} & \ldots & {{\bi x}_{k - 1}^T} \cr}} \right]^T $ and ${\bi w}_k = \left[ {\matrix{ {{\bi w}_{k1}} \cr {{\bi w}_{k2}} \cr}} \right]$, then Equation (26) is reformulated as:

(27)$${\bi x}_k = {\bi \Phi} _{k,\left( {k - n:k - 1} \right)} {\bi x}_{\left( {k - n:k - 1} \right)} + {\bi w}_k $$

where Φ_{k,(k − n:k − 1)} denotes the state transition matrix that is dependent on the time window length n and the time interval δt.

Once the window-length n and δt are given, Φ_{k,(k − n:k − 1)} will be a constant matrix defined as:

(28)$${\bi \Phi} _{k,\left( {k - n:k - 1} \right)} = \left[ {\matrix{ {{\bi \Phi} _{11}} & {{\bi \Phi} _{12}} \cr {{\bi \Phi} _{21}} & {{\bi \Phi} _{22}} \cr}} \right] = \left[ {\matrix{ {{\bi \psi} _{11}} & {{\bi \psi} _{12}} \cr {{\bi \psi} _{21}} & {{\bi \psi} _{22}} \cr}} \right] \otimes {\bi I}_{3 \times 3} $$

where ψ₁₁, ψ₁₂, ψ₂₁ and ψ₂₂ are all (n×n) coefficient matrices.

Their elements for the window length n from 1 to 5 are shown in Table 3.

Table 3. The coefficient matrix Ψ for the different window lengths (n = 1 to 5)

5. EXPERIMENT AND ANALYSIS

The experiment data were collected by using two Topcon HiPer-Pro GPS receivers with a sampling interval of 1 second. Although the maximal distance between the reference and rover stations is shorter than 3 km (see Figure 2), the various kinematic states are intentionally experienced in the whole vehicular GPS test. The estimated three-dimensional velocities of the vehicle by using the Doppler measurements are shown in Figure 3. The double differenced observation model is used in the experiment and the systematic errors are basically ignored due to their highly spatial correlation in such a short distance between reference and rover stations. The C/A code and Doppler D1 measurements were utilized in the whole test. The results computed by use of the dual-frequency carrier phase measurements (resolved by Ashtech solution GPS software 2.6) are used as the reference values to evaluate the accuracies of the developed models.

Figure 2. The horizontal (East – North) trajectory of vehicle.

Figure 3. Velocities of X, Y and Z components.

In the following section, we will investigate the accuracies of our WRA in the cases of the different window lengths and the different sampling intervals. In our computation strategy, the weight matrix of double differenced GPS observations was derived from the C/N₀ weighted matrix of undifferenced observations by using the law of error propagation (Kuusniemi et al., 2004). The researches have indicated that the stochastic model of w_k is also an important issue to affect the solution. Therefore, to obtain the optimal solution, the reasonable stochastic model should be determined in real-time (Wang 2000; Moore and Wang 2003, Hewitson and Wang 2007; Geng and Wang 2008; Yang and Gao 2005). In this paper, without affecting the validity and efficiency of our WRA, for simplification, the covariance matrix Σ _w was selected as a (6 × 6) matrix with spectral density of 0·2m²s⁻³ as follows (Schwarz et al., 1989):

(29)$${\bi \Sigma} _{{\bi w}_k} = \left[ {\matrix{ {\displaystyle{1 \over 3}\delta t^3} & {\displaystyle{1 \over 2}\delta t^2} \cr {\displaystyle{1 \over 2}\delta t^2} & {\delta t} \cr}} \right] \otimes {\bi Q}$$

where:

Q is a 3 × 3 diagonal matrix with the elements of spectral density for velocities.

δt denotes the sampling time interval.

The initial position and velocity variances were cautiously chosen as 1 m² and 0·1 m²/s² respectively. The transition matrix of the dynamic model established in Section 4 was used.

For the sampling rate of 1 s that is referred to as the high sampling rate, the differences of Least Squares (LS) solutions and WRA (n = 1) with/without Doppler measurements from reference values are shown in Figure 4. The accuracy statistics of LS and WRA (n = 1) solutions are presented in Table 4. It is easy to theoretically prove that the solution of WRA with n = 1 is exactly equivalent to the KF solution with CV dynamic model. Hence, herein they are not particularly compared with each other. Figure 5 shows the number of observed satellites (blue dashed line) and the Position Dilution Of Precision (PDOP) (red solid line). It can be seen from the Table 4 and Figure 4∼5 that:

1. (i) the poor solutions are generally assigned to the larger PDOP values, especially during the observation epochs from 2000 and 3000, which is mainly caused by the rapidly decreasing of observed satellites.

2. (ii) LS outperforms WRA (n = 1) without Doppler, because the kinematic model of WRA (n = 1) cannot accurately describe the moving object's manoeuvring motions, for instance, turning, accelerating and braking.

3. (iii) Aided by the kinematic model of the vehicle, compared with LS, WRA can improve the positioning accuracies in the case of the number of observed satellites sharply decreasing.

4. (iv) Since the Doppler measurements provide the accurate velocity quantities, WRA with Doppler is obviously better than that without Doppler, especially in poor observation conditions.

Figure 4. Differences of LS and WRA (n = 1) solutions from the reference values for sampling interval of 1 s (top: Differences of LS from the reference values; bottom: Differences of WRA without Doppler implementation).

Figure 5. The number of observed satellites and PDOP during kinematic GNSS navigation.

Table 4. RMSE (m) of WRA with and without Doppler implementations (n = 1) (1 s)

We also calculated the differences of the Doppler-aided WRA solutions based on the six-state model from reference values, which are illustrated in Figure 6. Moreover, the efficiency of constructed dynamic models of WRA is examined in the window length n from 1 to 5. The Root-Mean-Squared Error (RMSE) of three position components and position are computed by Equations (30) and (31) respectively:

(30)$${\rm RMS}\left( X \right) = \sqrt {\displaystyle{{\sum\nolimits_{i = 1}^q {{\rm\Delta} X_i^2}} \over q}} \quad {\rm RMS}\left( Y \right) = \sqrt {\displaystyle{{\sum\nolimits_{i = 1}^q {{\rm\Delta} Y_i^2}} \over q}} \quad {\rm RMS}\left( Z \right) = \sqrt {\displaystyle{{\sum\nolimits_{i = 1}^q {{\rm\Delta} Z_i^2}} \over q}} $$

(31)$${\rm RMS}\left( {POS} \right) = \sqrt {{\rm RMS}\left( X \right)^2 + {\rm RMS}\left( Y \right)^2 + {\rm RMS}\left( Z \right)^2} $$

where:

ΔX, ΔY and ΔZ represent the three position difference components between our WRA solutions and the referenced ones.

RMS(·) denotes the RMSE of position and the symbol POS denotes the point position.

Here q is the number of all observed epochs.

The computed RMSE statistics are shown in Table 5. Based on experimental results so far in Table 5 and Figure 6, it is concluded that:

1. (i) In the case of n = 2, the WRA solution is obviously superior to the solution of n = 1.

2. (ii) With the window length n increasing (n > 2), the accuracies decrease but the extent is not so obvious.

3. (iii) For those vehicles with high sampling rate (δt < 1 s) or manoeuvrability, the window length should not exceed 5 epochs. It is worth pointing out that the dynamic model in WRA of n = 2 actually works as the CA model and is also compatible with the CV model.

Figure 6. The differences of the WRA solutions from the reference values for sampling interval of 1 s.

Table 5. RMSE (m) of WRA with different window lengths (n = 1 to 5) in GPS navigation (1 s)

However, those dynamic models with larger windows (n > 2) do not achieve the expected solutions for two reasons:

1. (i) the dynamic characteristics are over modelled when n > 2.

2. (ii) The errors of NFDE and DI combined dynamic model sharply augmented with the window length n increasing, therefore the window length should be carefully chosen according to different application scenarios.

We further tested the performances of the Doppler-aided WRA in the lower sampling rate. The data is thinned in the sampling interval of 3 s and the computation strategy is the same as in the former experiment. The differences between the WRA solutions and the reference values are illustrated in Figure 7. The RMSE statistics computed in terms of Equations (30) and (31) are shown in Table 6. It is observed that:

1. (i) The RMSE statistics of WRA with different window lengths in Table 6 are generally greater than those in Table 5.

2. (ii) Some simple assumptions on the dynamic model of moving objects e.g. vehicle, CV (n = 1) and CA (n = 2), are not appropriate any more.

3. (iii) As the window length n increases, our WRA performs better and its accuracy is improved.

4. (iv) Compared with CV and CA dynamic models, those higher-order dynamic models (n > 2) accurately describe the characteristics of vehicle's motion and thus improve the navigation solutions in lower sampling rate scenario.

Figure 7. The differences of the WRA solutions from the reference values for sampling interval of 3 s.

Table 6. RMSE (m) of WRA with different window lengths (n = 1 to 5) in GPS navigation (3 s)

6. CONCLUSIONS

In summary, the implementation of Window-Recursive Approach (WRA) based on pseudorange and Doppler measurements has been developed for GNSS navigation. A state transition matrix was constructed based on the Newton's Forward Difference Extrapolation (NFDE) and Definite Integral (DI) methods for efficient computation. Real vehicular GPS experiments were carried out to evaluate the specified implementation of WRA in two sampling data scenarios. The results show that:

1. (i) in the high sampling rate, the low-order models can describe the motion states adequately and achieve the optimal results in the window length less than 5 epochs. Conversely, in the low sampling rate, the low order models cannot specify the motion states well enough. Therefore, high-order models with longer window length are necessary to achieve better performances.

2. (ii) WRA is flexible and efficient in Global Navigation Satellite Systems (GNSS) navigation computation. It reduces to a Kalman filter if window length n is set to 1.

3. (iii) The varying velocity and acceleration of a vehicle's motions can be automatically specified in the high-order models of WRA, therefore it does not need to construct the velocity and acceleration dynamic models for predicting the Position and Velocity (PV) state epoch by epoch.

In the future, we will further investigate the issues of WRA potential implementation fields, including optimal window length determination, nonlinear transformation problem and quality control framework etc. In addition, our WRA for GNSS kinematic navigation provides highly precise position and velocity information, which can also be implemented to accurately estimate Inertial Navigation System (INS) errors in GNSS/INS navigation or to detect faults in integrated navigation systems as a trajectory constraint condition.