2.1. Pseudorange measurement model

If we consider X_sat and X_recv to be the position vectors of the satellite and the receiver respectively such that in terms of the axial components x, y and z:

(1)$${\bi X}_{sat}=( x_{sat}\comma \; y_{sat}\comma \; z_{sat})$$

(2)$${\bi X}_{recv}=( x_{recv}\comma \; y_{recv}\comma \; z_{recv}) $$

then the geometric distance between the satellite and the receiver, ${}_{recv}^{sat} \rho $, is given by the l ₂-norm (Horn and Johnson, Reference Horn and Johnson2013) of the vector obtained by subtracting X_recv from X_sat:

(3)$${}_{recv}^{sat} \rho =\Vert {\bi X}_{sat}-{\bi X}_{recv}\Vert _{2} $$

Using the geometric distance ${}_{recv}^{sat} \rho $ obtained from Equation (3) and taking into account the receiver clock bias c _recvδt and satellite clock bias c ^satδt where c is the speed of light and δ t is one half of the time difference between sending and receiving the signal between the satellite and the receiver, along with the biases and delays associated with code range observables namely the receiver hardware bias ${}_{recv} \Delta_{P_{f}}$, satellite hardware bias ${}^{sat} \Delta_{P_{f}}$, ionospheric delay ${}_{recv}^{sat} I_{f}$, tropospheric delay ${}_{recv}^{sat} T_{f}$, multipath delay ${}_{recv}^{sat} m_{P_{f}}$ and other unmodelled errors $\varepsilon _{P_{f}}$, the Pseudorange measurement using code range observables, ${}_{recv}^{sat} P_{f}$, at a given frequency f can be modelled as:

(4)$${}_{recv}^{sat} P_{f}={}_{recv}^{sat} \rho +c( {}_{recv} {\delta t}-{}^{sat} {\delta t} ) +{}_{recv} \Delta_{P_{f}}-{}^{sat} \Delta_{P_{f}}+{}_{recv}^{sat} I_{f}+{}_{recv}^{sat} T+{}_{recv}^{sat} m_{P_{f}}+\varepsilon_{P_{f}} $$

The Pseudorange measurement at a given carrier frequency f using carrier phase observables ${}_{recv}^{sat} \Phi_{f}$ can also be modelled similarly to Equation (4) using equivalent bias and delay terms for carrier phase. However, in the case of the carrier phase, the ambiguous range also has to be accounted for, which is the ${}_{recv}^{sat} {\rm N}\lambda_{f}$ term in Equation (5), where N is the integer ambiguity and λ is the carrier wavelength.

(5)$$\eqalign{{}_{recv}^{sat} \Phi_{f}&={}_{recv}^{sat} \rho +c ( {}_{recv} {\delta t}-{}^{sat} {\delta t} ) +{}_{recv} \Delta _{\Phi_{f}}-{}^{sat} \Delta_{\Phi_{f}}+{}_{recv}^{sat} {{\rm N}\lambda}_{f} -{}_{recv}^{sat} I_{f} \cr & \quad +{}_{recv}^{sat} T+{}_{recv}^{sat} m_{\Phi_{f}}+\varepsilon_{\Phi_{f}}} $$

In the case of modelling the carrier phase measurement ${}_{recv}^{sat} \Phi_{f}$ as shown in Equation (5), the integer ambiguity N is a constant number as long as the carrier tracking loop maintains lock. Any break in tracking could change the value of N, which means that there is an increased complexity in assigning a reliable value for N in order to calculate the ambiguous range associated with carrier phase measurements (Mishra and Enge, Reference Mishra and Enge2004). Therefore, in order to keep the modelling process simple for the work presented in this paper, the Pseudorange measurement has been modelled using code range measurement, that is, Equation (4).

Equation (4) is further simplified by grouping together satellite-only, receiver-only and combined satellite-receiver errors so that the pseudorange measurement equation now becomes:

(6)$${}_{recv}^{sat} P^{\ast }={}_{recv}^{sat} \rho +{}_{recv} b-{}^{sat} b+{}_{recv}^{sat} b+\epsilon_{P^{\ast}} $$

These satellite-only, receiver-only and combined satellite-receiver errors are classified in Table 1. The measurement is then linearized about an approximate guess position for each satellite such that Equation (6) becomes:

Table 1.

Error classification (Anand, Reference Anand2018)

(7)$${}_{recv}^{sat} {\delta P}^{\ast }={}_{recv}^{sat} {\delta \rho }+{}_{recv} {\delta b}-{}^{sat} {\delta b}+{}_{recv}^{sat} {\delta b}+\epsilon_{P^{\ast }} $$

where δ represents the change in the value of each parameter with respect to that of the previous epoch. Further, the term ${}_{recv}^{sat} \delta \rho$ in Equation (7) can be modified using the Newton Raphson Method (Mishra and Enge, Reference Mishra and Enge2004) as:

(8)$${}_{recv}^{sat} {\delta \rho }=\left[\displaystyle{{x_{sat}-x_{recv}}\over{{}_{{\bi recv}}^{{\bi sat}} {\bi \rho}}}\comma \; \displaystyle{{y_{sat}-y_{recv}}\over{{}_{{\bi recv}}^{{\bi sat}} {\bi \rho}}}\comma \; \displaystyle{{z_{sat}-z_{recv}}\over{{}_{{\bi recv}}^{{\bi sat}} {\bi \rho}}} \right]^{{\bf T}}{\delta {\bi X}}_{sat} $$

In Equation (8),

(9)$$\delta {\bi X}_{sat}=( {\delta x}_{sat}\comma \; {\delta y}_{sat}\comma \; {\delta z}_{sat} ) $$

For the sake of simplicity let:

(10)$$\left[\displaystyle{{x_{sat}-x_{recv}}\over{{}_{{\bi recv}}^{{\bi sat}} {\bi \rho}}}\comma \; \displaystyle{{y_{sat}-y_{recv}}\over{{}_{{\bi recv}}^{{\bi sat}} {\bi \rho}}}\comma \; \displaystyle{{z_{sat}-z_{recv}}\over{{}_{{\bi recv}}^{{\bi sat}} {\bi \rho}}} \right]^{{\bf T}}=( -{{\bf 1}}_{{\bi recv}}^{{\bi sat}} ) $$

So that Equation (7) can be re-written as:

(11)$${}_{recv}^{sat} {\delta P}^{\ast }=( -{{\bf 1}}_{{\bi recv}}^{{\bi sat}} ) {\delta {\bi X}}_{sat}+{}_{recv} {\delta b}-{}^{sat} {\delta b}+{}_{recv}^{sat} {\delta b}+\epsilon_{P^{\ast }} $$

For n satellites and m receivers, we obtain mn number of pseudorange measurements. These measurements are then linearized as shown in Equation (11) to build a system of equations as shown in Equation (12):

(12)$${\bf A}_{mn\times 3n}{\bf u}^{\ast} +{\bf I}_{mn}{\bf u}_{{\rm los}}+{\bf B}_{mn\times ( m+n) }{\bf u}_{{\rm elim}}={\rm y} $$

In this system of linearized equations, y is a matrix of size mn × 1 and it contains all the values of ${}_{recv}^{sat} {\delta P}^{\ast }$ for n satellites and m receivers. I is an identity matrix of size mn × mn, u is the vector containing the position, biases and the delays split into: u* which contains change in position vectors over an iteration for each satellite:

(13)$${\bf u}^{\ast }=[ {}^{1} ( \delta x_{sat}\ \delta y_{sat}\ \delta z_{sat}) \comma \; \ldots^{n} ( \delta x_{sat}\ \delta y_{sat}\ \delta z_{sat}) ] _{3n\times 1}^{\rm T} $$

u_los contains the ionospheric, tropospheric and multipath delays which can be modelled from measurements and removed:

(14)$${\bf u}_{{\rm los}}={[ {}_{1}^{1} {\delta b}\comma \; \ldots {}_{m}^{n} {\delta b}] }_{mn\times 1}^{\rm T} $$

and u_elim contains the satellite-only and receiver-only biases:

(15)$${\bf u}_{{\rm elim}}={[ {}_{1} {\delta b}\comma \; \ldots\ {}_{m}{\delta b}\comma \; {}^{1}{\delta b}\comma \; \ldots {}^{n}{\delta b}] }_{( m+n) \times 1}^{\rm T} $$

A is the Jacobian of the pseudorange measurements with respect to the position vector as shown in Equation (16):

(16)$${\bf A}=\left[\matrix{( {\bf -1}_{{\bf 1}}^{{\bf 1}})^{{\bf T}} & {\open O}_{{\bf 3}}^{{\bf T}} & \cdots & {\open O}_{{\bf 3}}^{{\bf T}} \cr \vdots & \vdots & \, \, & \vdots \cr ( {\bf -1}_{{\bi m}}^{{\bf 1}})^{{\bf T}} & {\open O}_{{\bf 3}}^{{\bf T}} & \cdots & {\open O}_{{\bf 3}}^{{\bf T}} \cr {\open O}_{{\bf 3}}^{{\bf T}} & ( {{\bf -1}}_{{\bf 1}}^{{\bf 2}})^{{\bf T}} & \cdots & {\open O}_{{\bf 3}}^{{\bf T}} \cr \vdots & \vdots & \, \, & \vdots \cr {\open O}_{{\bf 3}}^{{\bf T}} & ( {\bf -1}_{{\bi m}}^{{\bf 2}})^{{\bf T}} & \, \, & {\open O}_{{\bf 3}}^{{\bf T}} \cr \vdots & \vdots & \, \, & \vdots \cr {\open O}_{{\bf 3}}^{{\bf T}} & {\open O}_{{\bf 3}}^{{\bf T}} & \cdots & ( {{\bf -1}}_{{\bf 1}}^{{bi n}})^{{\bf T}} \cr \vdots & \vdots & \, \, & \vdots \cr {\open O}_{{\bf 3}}^{{\bf T}} & {{\open O}}_{{\bf 3}}^{{\bf T}} & \cdots & ( {{\bf -1}}_{{\bi m}}^{{bi n}})^{{\bf T}}}\right]_{mn\times 3n} $$

In Equation (16), the terms $( {{\bf -1}}_{{\bi m}}^{{bi n}}) $ are as defined in Equation (10) while 𝟘₃ is a 3 × 1 matrix of zeros.

B is the Jacobian of the pseudorange measurements with respect to the satellite-only and receiver-only biases:

(17)$${\bf B}=\left[\matrix{{\bf I}_{m} & -{\open I}_{m} & {\open O}_{m} & \ldots & {\open O}_{m} \cr {\bf I}_{m} & {\open O}_{m} & -{\open I}_{m} & \ldots & {\open O}_{m} \cr \vdots & \vdots & \vdots & \, \, & \vdots \cr {\bf I}_{m} & {\open O}_{m} & {\open O}_{m} & \ldots & -{\open I}_{m}}\right]_{mn\times ( m+n) } $$

In Equation (17), I is an identity matrix of size m × m, 𝟙_m is a m × 1 matrix of ones and 𝟘_m is a m × 1 matrix of zeros.

2.2. Projection method

The Single Difference method eliminates the satellite-only bias terms from the measurement function, while the Double Difference (DD) method goes one step further and eliminates the receiver-only bias terms as well, thus mathematically shifting focus towards the position estimation problem without having to account for satellite and receiver-only biases (Mishra and Enge, Reference Mishra and Enge2004). The downside to this is that double differencing reduces the number of available observations while causing higher error correlation and noise.

De Jonge (Reference De Jonge1998) shows a method of undifference that helps to partly counter the above problem using a projection matrix. The projection method allows for the number of available observations to remain unaffected while reducing the degree of error correlation with increasing satellite and receiver combinations although this results in a marginal increase in noise (Anand, Reference Anand2018).

For the current system of n satellites and m receivers, there exists a K1 and K2 such that,

(18)$$K1=\left[\matrix{( m-1) ( n-1) & 1-n & \ldots & 1-n \cr 1-n & ( m-1) ( n-1) & \ldots & 1-n \cr \vdots & \vdots & \, \, & \vdots \cr 1-n & 1-n & \ldots & ( m-1) ( n-1) \cr }\right]_{m\times m} $$

(19)$$K2=\left[\matrix{1-m & 1 & \cdots & 1 \cr 1 & 1-m & \cdots & 1 \cr \vdots & \vdots & \, \, & \vdots \cr 1 & 1 & \cdots & 1-m}\right]_{m\times m} $$

(20)$${\bf K}=\displaystyle{{1}\over{mn}}\left[\matrix{K1 & K2 & \cdots & K2 \cr K2 & K1 & \cdots & K2 \cr \vdots & \vdots & \, \, & \vdots \cr K2 & K2 & \cdots & K1}\right]_{mn\times mn} $$

where K is the projection matrix for the Jacobian B. Therefore, we arrive at a helpful result: KB = 0_mn×(m+n).

The projection method thus allows for the nullification of the satellite-only and receiver-only bias vectors. The updated linearized measurement equation now becomes:

(21)$${\bf KAu}^{\ast}={\bf Ky}\Rightarrow {\bf A}^{\prime}{\bf u}^{\ast}={\bf y}^{\prime} $$

The measurement covariance is updated as:

(22)$${\brSigma}_{{\bf y}^{\prime}}=\sigma_{P} \ast {\bf KK}^{\rm T}=\sigma_{P}\ast {\bf K} $$

since the property of the projection matrix is such that K^T = K and K² = K (Strang, Reference Strang2010).

u* is then estimated by calculating the Moore-Penrose inverse (Penrose, Reference Penrose1955) of A^′:

(23)$${\bf u}^{\ast}=( {\bf A}^{\prime}) ^{+}{\bf y}^{\prime} $$

(24)$${\brSigma}_{{\bf u}}^{\ast}=( {\bf A}^{\prime}{\brSigma}_{{\bf y}^{\prime}} {\bf A}^{\prime{\rm T}}) ^{-1} $$

2.3. Kalman filter formulation

In this section, a Modified Kalman Filter (MKF) is formulated in a similar fashion to the familiar Conventional Extended Kalman Filter (CKF) algorithm and following the projection-based measurement function formulated in the previous section.

Table 2 summarises the algorithms for both the CKF and the MKF and provides a comparison of the formulation of both filters.

Table 2.

Summary of CKF and MKF algorithms

In both the CKF and the MKF, X represents the satellite's State Vector (SV):

(25)$${\bf X}=[ x\ y\ z\ \dot{x}\ \dot{y}\ \dot{z}] _{6\times 1}^{\rm T} $$

Δ X is the update in state vector added to the SV from the previous epoch, X_i−1:

(26)$${\brDelta}{\bf X}=[ \delta x\ \delta y\ \delta z\ \dot{\delta x}\ \dot{\delta y}\ \dot{\delta z}] _{6\times 1}^{\rm T} $$

In Table 2, Φ is the dependence of the SV on the initial state, that is, X at t = 0. However, unlike a least squares estimation algorithm that uses all the known measurements to estimate position at every epoch, the Kalman Filter holds an advantage in being an online estimation algorithm that follows an incremental fusing of estimates (Pei et al., Reference Pei, Biswas, Fussell and Pingali2017). Therefore, the state transition matrix in this case is the dependence of the SV on the state at the previous epoch that is:

(27)$${\brPhi}_i=\displaystyle{{\partial {\bi X}_i} \over {\partial {\bi X}_{i-1}}} $$

P is the covariance of the position estimate as formulated in the previous section and K is the Kalman gain matrix. In both algorithms shown in Table 2, a minus sign in the superscript indicates an update in the model function while a plus sign indicates an update in the measurement function. As is the case in any Kalman Filter algorithm, the measurement update from the previous epoch serves as the a priori estimate for the model update in the succeeding epoch.

The difference between the algorithms is that the MKF utilises the projection-based position estimates from the measurement function while the CKF includes the satellite and receiver biases so that the innovation in the CKF is characterised by:

(28)$${\bf y}=[ {}_{{\bf 1}}^{{\bi sat}} {\delta P}\comma \; {}_{{\bf 2}}^{{\bi sat}} {\delta P}\comma \; \ldots {}_{{\bi m}-{\bf 1}}^{{\bi sat}} {\delta P}\comma \; {}_{{\bf m}}^{{\bi sat}} {\delta P}]_{{\bi m}\times {\bf 1}}^{\rm T} $$

where P is the pseudorange measurement between the satellite and the receiver. The Modified Kalman Filter on the other hand utilises Equation (29) which is the output from the projection method:

(29)$${\bf u}=[ \delta x\, \delta y\, \delta z] _{3\times 1}^{\rm T} $$