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 Reference Brown and Hwang1997; Simon Reference Simon2006). 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 Reference Julier and Uhlmann1997; Julier et al., Reference Julier, Uhlmann and Durrant-Whyte2000; Jwo and Lai Reference Jwo and Lai2009) 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 Reference Jwo and Lai2008).

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., (Reference Zhou, Shen and Li2010). 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 Reference Han and Wang2012). 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., Reference Yang, Cui and Gao2004).

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., Reference Zhou, Shen and Li2010). 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., (Reference Zhou, Li, Rizos and Shen2009, Reference Zhou, Shen and Li2010), 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., Reference Zhou, Shen and Li2010):

(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., Reference Kuusniemi, Lachapelle and Takala2004). 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 Reference Wang2000; Moore and Wang Reference Moore and Wang2003, Hewitson and Wang Reference Hewitson and Wang2007; Geng and Wang Reference Geng and Wang2008; Yang and Gao Reference Yang and Gao2005). 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., Reference Schwarz, Cannon and Wong1989):

(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)