2.1 GNSS measurement model

A GNSS receiver measures the distance between the satellite and the receiver's antenna by multiplying the speed of light by the propagation time. Any clock-related errors and propagation errors caused by the ionosphere and troposphere are included in the measured distance, and the measured range differs from the true range. Thus, the basic pseudorange measurement is written as (Groves, Reference Groves2008):

(1)\begin{equation}{\rho _i}= \sqrt {{{(x_i^s-{x^r})}^2} + {{(y_i^s-{-}{y^r})}^2} + {{(z_i^s-{z^r})}^2}} + {I_{\textrm{trop}}} + {I_{\textrm{iono}}} + c\, \cdot\, d{t^r} - c \,\cdot\, d{t^s} + \partial\end{equation}

where ${\rho _i}$ denotes the measured pseudorange; $(x_i^s,y_i^s,z_i^s)$ and $({x^r},{y^r},{z^r})$ denote the position of the ${i^{th}}$ satellite and the receiver's antenna in the Earth-centred and Earth-fixed (ECEF) frame, respectively; ${I_{\textrm{iono}}}$ and ${I_{\textrm{trop}}}$ are the ionospheric delay and tropospheric delay, respectively; $d{t^r}$and $d{t^s}$ denote the clock bias of the receiver and satellite, respectively; c represents the speed of light; and $\partial$ represents the remaining error, including the noise.

The pseudorange rate refers to the velocity obtained by differentiating the pseudo-distance in the measurement time interval and could transform from the measured Doppler shift. The basic pseudorange rate measurement could be expressed as (Groves, Reference Groves2008):

(2)\begin{equation}{\dot{\rho }_i}= {\bf D} \,\cdot\, [{{\bf v}_i^s - {{\bf v}^r}} ]+ {\dot{I}_{\textrm{trop}}} + {\dot{I}_{\textrm{iono}}} + c \,\cdot\, \dot{d}{t^r} - c \,\cdot\, \dot{d}{t^s} + \dot{\partial }\end{equation}

where ${\dot{\rho }_i}$ denotes the measured pseudorange rate; ${\bf D}$ denotes the line-of-sight vector from the ${i^{th}}$satellite to the GNSS receiver's antenna; ${\bf v}_i^s$ and ${{\bf v}^r}$denote the velocity vector of the ${i^{th}}$ satellite and the receiver, respectively; ${\dot{I}_{\textrm{trop}}}$ and ${\dot{I}_{\textrm{iono}}}$ are tropospheric rate correction and ionospheric rate correction, respectively; $\dot{d}{t^r}$ and $\dot{d}{t^s}$ denote the clock drift of the receiver and satellite, respectively; and $\dot{\partial }$ represents the remaining errors, including noise.

2.2 INS Mechanisation

An INS includes a navigation computer and inertial measurement unit (IMU), which is a self-contained navigation system. The accelerometers and gyroscopes from the IMU are used to measure the user acceleration and angular velocity in body frame (Right-Front-Up, RFU). The INS mechanisation is implemented in the navigation frame (East-North-Up, ENU) to update position, velocity and attitude. The navigation solution updating equations are given as (Groves, Reference Groves2008):

(3)\begin{gather}{\dot{{\bf r}}_{\textrm{IN}{\textrm{S}_n}}} = {{\bf v}_n} - \boldsymbol{\omega}_{en}^n{{\bf r}_{\textrm{IN}{\textrm{S}_n}}}\end{gather}

(4)\begin{gather}{\dot{{\bf v}}_{\textrm{IN}{\textrm{S}_n}}} = {{\bf f}_n} + {{\bf g}_n} - (\boldsymbol{\omega}_{en}^n + 2\boldsymbol{\omega}_{ie}^n) \times {{\bf v}_{\textrm{INS}}}_{_n}\end{gather}

(5)\begin{gather}\dot{{\bf q}}_b^n = \frac{1}{2}{\bf q}_b^n\boldsymbol{\omega}_{nbq}^b\end{gather}

where the subscripts/superscripts ‘e’, ‘n’ and ‘i’ represent the ECEF frame, navigation frame and earth-centred inertial frame, respectively; ${{\bf r}_{\textrm{IN}{\textrm{S}_n}}}$ and $\dot{{\bf r}}_{\textrm{IN}{\textrm{S}_n}}$ represent the position vector and corresponding position rate vector, respectively;${{\bf v}_{\textrm{INS}}}_{n}$ and ${\dot{{\bf v}}_{\textrm{IN}{\textrm{S}_n}}}$ represent the velocity vector and corresponding velocity rate vector respectively; ${{\bf g}_n}$ refers to the acceleration vector; ${{\bf f}_n}$ represents the specific force vector; $\boldsymbol{\omega}_{ie}^n$ and $\boldsymbol{\omega}_{en}^n$ denote the angular velocity vector with respect to the inertial frame and the earth fixed frame respectively; ${\bf q}_b^n$ and $\dot{{\bf q}}_b^n$ are the attitude quaternions and its attitude quaternion rate; and $\boldsymbol{\omega}_{nbq}^b$ represents the skew-symmetric matrix of the angular velocity.

In practical applications, the lever arm between the GNSS receiver's antenna to the IMU must be considered due to the different install position of the INS device and the GNSS receiver's antenna. In this paper, we chose to correct the centre position of the INS via the INS-derived position and attitude solution to the GNSS antenna position. The INS-derived position at the GNSS antenna centre is calculated as:

(6)\begin{equation}{{\bf r}_{\textrm{INS}\_\textrm{GNS}{\textrm{S}_n}}} = {{\bf r}_{\textrm{IN}{\textrm{S}_n}}} + {\bf C}_b^n{{\bf L}_b}\end{equation}

where ${{\bf r}_{\textrm{IN}{\textrm{S}_n}}}$ is the INS-derived position at the IMU centre; ${{\bf r}_{\textrm{INS}\_\textrm{GNS}{\textrm{S}_n}}}$ is the corrected INS position at the GNSS receiver's antenna; ${{\bf L}_b}$ is the lever arm vector; and ${\bf C}_b^n$ represents the attitude transformation matrix.

2.3 Error-state based TC GNSS/INS integration

In the TC GNSS/INS integration, the raw GNSS pseudorange rate and pseudorange are fused with INS-derived pseudorange rate and pseudorange. The system filter uses input measurements (pseudorange rate difference and pseudorange difference) to estimate error states to correct the corresponding INS-derived absolute states to obtain corrected absolute states.

The system state vector ${\bf x}(t)$of the TC integration method contains seventeen states and is given as:

(7)\begin{equation}{\bf x}(t) = {[\delta {r_E},\delta {r_N},\delta {r_U},\delta {v_E},\delta {v_N},\delta {v_U},{\varphi _E},{\varphi _N},{\varphi _U},\textrm{ }{a_x},{a_y},{a_z},{\eta _x},{\eta _y},{\eta _z},\delta d{t^r},\delta \dot{d}{t^r}]^T}\end{equation}

where $(\delta {r_E},\delta {r_N},\delta {r_U})$ and $(\delta {v_E},\delta {v_N},\delta {v_U})$ represent position errors and velocity errors in the navigation frame; $({\varphi _E},{\varphi _N},{\varphi _U})$ represents attitude error in the body frame; $({a_x},{a_y},{a_z})$ and $({\eta _x},{\eta _y},{\eta _z})$ represent accelerometer bias error and gyroscope drift error; $\delta d{t^r}$ and $\delta \dot{d}{t^r}$ represent clock bias error and clock drift error of the receiver.

The system measurement consists of pseudorange and pseudorange rate difference, and the measurement vector can be written as:

(8)\begin{equation} {\bf y}(t)= {[{{\boldsymbol{\rho}_{\textrm{GNSS}}} - {\boldsymbol{\rho}_{\textrm{INS}}}\textrm{ }{{\dot{\boldsymbol{\rho}}}_{\textrm{GNSS}}} - {{\dot{\boldsymbol{\rho}}}_{\textrm{INS}}}} ]^T}\end{equation}

where ${\boldsymbol{\rho}_{\textrm{GNSS}}}$ and ${\boldsymbol{\rho}_{\textrm{INS}}}$ are the GNSS raw pseudorange and INS computed pseudorange vectors, respectively; and ${\dot{\boldsymbol{\rho}}_{\textrm{GNSS}}}$ and ${\dot{\boldsymbol{\rho}}_{\textrm{INS}}}$ are the GNSS raw pseudorange rate and INS computed pseudorange rate vectors, respectively.

The widely used tightly coupled GNSS/INS is based on the EKF. The discrete time system dynamic model and measurement model are:

(9)\begin{gather}{{\bf x}_k} = {\boldsymbol{\Phi}_k}{{\bf x}_{k - 1}} + {\boldsymbol{\varsigma}_k}\end{gather}

(10)\begin{gather}{{\bf y}_k} = {{\bf H}_k}{{\bf x}_k} + {{\bf o}_k}\end{gather}

where ${{\bf x}_k}$ and ${{\bf x}_{k - 1}}$ represent the system state vector at k and $k - 1$ epoch; ${\boldsymbol{\Phi}_k}$ and ${\boldsymbol{\varsigma}_k}$ represent the state transition matrix and system processing noise; ${{\bf y}_k}$ and ${{\bf H}_k}$ represent the measurement vector and measurement matrix; and ${{\bf o}_k}$ represents the measurement noise. Details of these vectors and matrixes are included in our previous research (Wei et al., Reference Wei, Dan, Baigen, Chris, Jian and Wei2019).

Then the standard Kalman filter is then applied to update time and measurement.

3.1 TC GNSS/INS with error-state UKF

The conventional EKF can estimate the system state in a tightly coupled GNSS/INS integration method, but it brings out truncation, resulting in reduced estimation accuracy and high computational efforts by expanding at an approximate point with first-order Taylor.

Differing from the EKF, the UKF adopts propagated sigma points to generate transformed sigma points, and then the transformed sigma points are weighted and further summed to estimate the expectation. This can overcome the problem of linearised truncation to provide third-order estimation accuracy, while the EKF can only result in first-order accuracy. Hence, the UKF is preferred in the research and the TC GNSS/INS integration with error-state UKF is designed to provide navigation solutions in GNSS open-sky areas. The error-state UKF-based tightly coupled GNSS/INS method consists of three steps: time update, unscented transformation and measurement update.

Due to the linearity of the dynamic model for GNSS/INS integration system in this research (Goshen-Meskin and Bar-Itzhack, Reference Goshen-Meskin and Bar-Itzhack1992), the time update of UKF can be simplified and conducted as that of standard Kalman filter:

(11)\begin{gather}\hat{{\bf x}}_k^ -= {\boldsymbol{\Phi}_{k - 1}}{\hat{{\bf x}}_{k - 1}}\end{gather}

(12)\begin{gather}{\bf P}_k^ -= {\boldsymbol{\Phi}_{k - 1}}{{\bf P}_{k - 1}}{\boldsymbol{\Phi}_{k - 1}}^T + {{\bf Q}_{k - 1}}\end{gather}

The key of UKF is using UT to generate the sigma points of an n-dimensional random variable, and associated weights. In this research, the one-step predicted error state $\hat{{\bf x}}_k^ -$ is regarded as the n-dimensional random variable, the corresponding mean and covariance are $\hat{{\bf x}}_k^ -$ and ${\bf P}_k^ -$. Then, the 2n + 1 sigma vector $\hat{\boldsymbol{\chi}}$ is calculated in detail as:

(13)\begin{gather}{\begin{array}{*{20}{c}} {\hat{\boldsymbol{\chi}}_k^ - } \end{array}_0} = \hat{{\bf x}}_k^ -\end{gather}

(14)\begin{gather}{\begin{array}{*{20}{c}} {\hat{\boldsymbol{\chi}}_k^ - } \end{array}_j} = \hat{{\bf x}}_k^ -+ {\left( {\sqrt {(n + \zeta ){\bf P}_k^ - } } \right)_j},j = 1, \cdots ,n\end{gather}

\begin{gather}{\begin{array}{*{20}{c}} {\hat{\boldsymbol{\chi}}_k^ - } \end{array}_{j + n}} = \hat{{\bf x}}_k^ - \textrm{ - }{\left( {\sqrt {(n + \zeta ){\bf P}_k^ - } } \right)_j},j = 1, \cdots ,n\nonumber\end{gather}

where ${\left( {\sqrt {(n + \zeta ){\bf P}_k^ - } } \right)_j}$ represents the j ^th column of the matrix square root; n denotes the dimension of the error state vector ${\bf x}(t)$ (it is 17 in this paper); $\zeta = {\alpha ^2}(n + \lambda ) - n$ is a scaling parameter determined by the dimension of the error state vector n and two tuning parameters $\alpha$ and $\lambda$; $\alpha$ determines the spread of the sigma points around $\hat{{\bf x}}_k^ -$ (set to $1e - 3$ in this paper); and $\lambda$ is usually set as 0 (Julier et al., Reference Julier, Uhlmann and Durrant-Whyte1995).

The associate weights ${\bf G}$ is given as:

(15)\begin{gather}{\bf G}_0^{\textrm{mean}} = \frac{\zeta }{{n + \zeta }}\end{gather}

(16)\begin{gather}{\bf G}_0^{{\mathop{\rm cov}} } = {\bf G}_0^{\textrm{mean}} + (1 - {\alpha ^2} + b)\end{gather}

\begin{gather}{\bf G}_j^{mean} = {\bf G}_j^{{\mathop{\rm cov}} } = \frac{1}{{2(n + \zeta )}},j = 1,\ldots ,2n\nonumber\end{gather}

where ${\bf G}_j^{\textrm{mean}}$and ${\bf G}_j^{{\mathop{\rm cov}} }$ are the mean weight and covariance weight associated with the j ^th point, respectively; and b is used to incorporate any a priori knowledge about the distribution of $\hat{{\bf x}}_k^ -$and the optimal value is 2 (Julier et al., Reference Julier, Uhlmann and Durrant-Whyte1995).

The nonlinear measurement model of error-state UKF is expressed as:

(17)\begin{equation}{\bf y}(t) = {\bf h}({\bf x}(t),{\bf o}(t))\end{equation}

specifically,

(18)\begin{equation}{\bf h}({\bf x}(t)) = [{{\bf h}_\rho }({\bf x}(t))\textrm{ }{{\bf h}_{\dot{\rho }}}({\bf x}(t))\textrm{]}\end{equation}

where ${{\bf h}_\rho }({\bf x}(t))$ and ${{\bf h}_{\dot{\rho }}}({\bf x}(t))$ are the pseudorange and pseudorange rate measurement functions in Equations (1) and (2).

Note that the pseudorange rate measurement model in Equation (2) is linear in theory. To maintain the simplicity and consistency of the system measurement, the original pseudorange rate measurement model performs 2n + 1 linear computations using 2n + 1 sigma points and is used as a nonlinear measurement model in this research.

Substituting the sigma points and associate weights of the one-step predicted error state into Equation (17), the prediction measurement points are created as:

(19)\begin{equation}{\begin{array}{*{20}{c}} {{\bf y}_k^ - } \end{array}_j} = {\bf h}[{\begin{array}{*{20}{c}} {\hat{\boldsymbol{\chi}}_k^ - } \end{array}_j}],\,j = 0,\ldots ,2n\end{equation}

The mean of the predicted measurement is computed as:

(20)\begin{equation}\hat{{\bf y}}_k^ -= \sum\limits_{j = 0}^{2n} {{\bf G}_j^{\textrm{mean}}{{\begin{array}{*{20}{c}} {{\bf y}_k^ - } \end{array}}_j}}\end{equation}

Then the measurement innovation covariance and the cross-covariance are calculated as:

(21)\begin{gather}\begin{array}{*{20}{c}} {{\bf P}_{zz}^-}_k \end{array} = \sum\limits_0^{2n} {{\bf G}_j^{{\mathop{\rm cov}} }[{{\begin{array}{*{20}{c}} {{\bf y}_k^ - } \end{array}}_j} - \hat{{\bf y}}_k^ - ] \,\cdot\, } {[{\begin{array}{*{20}{c}} {{\bf y}_k^ - } \end{array}_j} - \hat{{\bf y}}_k^ - ]^T} + {{\bf R}_k}\end{gather}

(22)\begin{gather}{{\bf P}_{xz}^-}_k= \sum\limits_0^{2n} {{\bf G}_j^{{\mathop{\rm cov}} }[{{\begin{array}{*{20}{c}} {\hat{\boldsymbol{\chi}}_k^ - } \end{array}}_j} - \hat{{\bf x}}_k^ - ] \,\cdot\, } [{\begin{array}{*{20}{c}} {{\bf y}_k^ - } \end{array}_j} - \hat{{\bf y}}_k^ - ]\end{gather}

The gain, state and corresponding state error covariance are updated as:

(23)\begin{gather}{{\bf K}_k}= {{\bf P}_{xz}^-}_k \begin{array}{*{20}{c}} {{\bf P}_{zz}^-}_k \end{array}\end{gather}

(24)\begin{gather}\hat{{\bf x}}_k^ += \hat{{\bf x}}_k^ -+ {{\bf K}_k}[{{\bf y}_k} - \hat{{\bf y}}_k^ - ]\end{gather}

\begin{gather}{\bf P}_k^ += {\bf P}_k^ - \textrm{ }{{\bf K}_k}\begin{array}{*{20}{c}} {{\bf P}_{zz}^-}_k \end{array}{\bf K}_k^T\nonumber\end{gather}

3.2 Advanced MM-augmented TC GNSS/INS with error-state UKF

In practice, the train usually operates in GNSS partly blocked areas, such as mountains, valleys, railway stations and so on. In such environments, the visible satellites are reduced and the GNSS geometry becomes worse, which results in poor GNSS measurement accuracy. As a result, the positioning performance of the TC GNSS/INS integration method also decreases. To maintain the accuracy and robustness of the integration system, an advanced MM-augmented TC GNSS/INS integration method with error-state UKF is proposed.

A typical MM application consists of two steps: search and projection, as shown in Figure 1. The first step is to search for the two consecutive key points closest to the original train location from all key points of the map. Due to the close distance between these two consecutive key points, the line segment between the two of them can be approximated as a straight line, called the nearest match interval. Then the original train location is projected perpendicularly on the nearest match interval to obtain a new point. The new point is called the map-matching corrected train location and output as the final train location.

Figure 1. Typical map-matching steps

Different from the traditional MM method with the shortest path search algorithm (Wei et al., Reference Wei, Sirui, Baigen, Jian and Wei2018), the proposed advanced MM-augmented TC GNSS/INS integration method uses a binary search algorithm to quickly search for the nearest interval with less comparison times to improve the operation efficiency. In addition, the original train location is derived from the one-step predicted location not the final estimated location after the measurement update. This altered original train location can avoid poor positioning results due to poor measurement quality caused by GNSS signal blockage. Most important of all, the MM-corrected location further fuses with INS as new measurement to improve estimation accuracy instead of output directly. This proposed method is able to improve positioning accuracy and generate smoother positioning solutions in GNSS partly blocked areas.

3.2.1 Binary search algorithm

In the binary search algorithm, the first step is to determine the upper and lower search index u and l for the original point. As the binary search algorithm begins at both the ends of the DTM, the initial value of the upper search index is set to 1 and the total number of track points. Then, the corresponding points are assumed as ${L_u}\textrm{ }({{M_u},{N_u}} )$ and ${L_l}\textrm{ }({{M_l},{N_l}} )$.

Define the original point as $K({M,N} )$ and the distances from point $K({M,N} )$ to the two current search points as ${d_u}$ and ${d_l}$. The corresponding distance calculation function is:

(25)\begin{equation}\left\{ {\begin{array}{*{20}{c}} {{d_u} = \sqrt {{{(M - {M_u})}^2} + {{(N - {N_u})}^2}} }\\ {{d_l} = \sqrt {{{(M - {M_l})}^2} + {{(N - {N_l})}^2}} } \end{array}} \right.\end{equation}

At the same time, the middle search index is updated by averaging the upper and lower search indexes as:

(26)\begin{equation}\textrm{mid} = \textrm{fix}[{(u + l)/2} ]\end{equation}

The upper and lower search indexes are constantly updated by comparing the values of ${d_u}$ and ${d_l}$. Assign the value of middle search index to the lower search index if ${d_u}$ is great than ${d_l}$. Otherwise, the value of middle search index is assigned to the upper search index. Then the nearest point ${L_{\textrm{mid}}}\textrm{ }({{M_{\textrm{mid}}},{N_{\textrm{mid}}}} )$ will be determined if the lower index u is equal or greater than the upper index l. The other nearest point (${L_{\textrm{mid} - 1}}$ or ${L_{\textrm{mid + }1}}$) is then located, which is the closest point to $K({M,N} )$ between the left and right points of the nearest point ${L_{\textrm{mid}}}$. Based on the two nearest points, the nearest match interval is determined.

3.2.2 Advanced MM measurement-augmented method

The accuracy of the MM depends on the accuracy of the DTM and the original point. The DTM is usually generated by a series of highly precise and consecutive track points. To ensure the MM accuracy, the distance between two adjacent points of the DTM is about 3 meters in this paper. For the traditional MM-aided TC GNSS/INS method, the original point is the TC integration position (final estimated position from UKF measurement step). However, the accuracy of the integration position is heavily degraded in GNSS partly blocked areas, and the traditional MM-aided method is unable to ensure the matching accuracy.

To solve the problem, this research considers MM before the degraded GNSS fusion step, and moves MM forward from the measurement update step to the time update step of the UKF. In the time update step, the one-predicted position can be calculated by one-step predicted position error and the INS-derived position in a short time. Note that the one-step predicted position error position is usually set to zero due to the closed-loop correction in the TC GNSS/INS integration method. Thus, the accuracy of the one-step predicted position mainly depends on the INS-derived position, which performs short-time accuracy. In this paper, the one-step predicted position outputs as an original point for MM and the corresponding matched position is generated via point-line MM. Due to the short-time accuracy of the matched position, it can be taken as a new measurement to fuse with INS to augment TC GNSS/INS instead of direct output in GNSS partly blocked areas. This matched position is defined as an advanced MM measurement.

The one-step predicted point is perpendicularly projected on the nearest match straight to obtain the MM point ${{\bf r}_{\textrm{MM}}}$. An additional position is updated by fusing n advanced MM position with the INS-derived position in the position domain. Hence, the addition measurement vector is calculated by advanced MM position and INS position difference and is written as:

(27)\begin{equation}{{\bf y}_{\textrm{MM}}}(t) = [{{\bf r}_{\textrm{MM}}} - {{\bf r}_{\textrm{INS}}}]\end{equation}

The additional measurement model can be described as:

(28)\begin{equation}{{\bf y}_{\textrm{MM}}}(t) = {{\bf H}_{\textrm{MM}}}(t){\bf x}(t) + {{\bf o}_{\textrm{MM}}}(t)\end{equation}

where ${{\bf o}_{\textrm{MM}}}(t)$ denotes measurement noise of the MM point and the corresponding covariance matrix is ${{\bf R}_{\textrm{MM}}}$; and ${{\bf H}_{\textrm{MM}}}(t)$ denotes the measurement matrix and is calculated as:

(29)\begin{equation}{{\bf H}_{\textrm{MM}}} = \left[ {\begin{array}{*{20}{c}} {{{\bf I}_{3 \times 3}}}& {{{\bf 0}_{3 \times 3}}}& {{{\bf 0}_{3 \times 3}}}& {{{\bf 0}_{3 \times 3}}}& {{{\bf 0}_{3 \times 3}}}& {{{\bf 0}_{2 \times 2}}} \end{array}} \right]\end{equation}

With the poor GNSS pseudorange rate and pseudorange in GNSS partly blocked areas, the accuracy of integration navigation solutions is decreased. The purpose of this research is to take advantage of available high-precision geographic information on DTM without adding additional sensors. Hence, the new MM position from additional position updates is applied to aid GNSS measurements and enhance the measurement model to maintain system positioning accuracy.

Compared to the standalone UKF-based tightly coupled GNSS/INS integration method, only the measurement vector and measurement model of the advanced MM-enhanced tightly coupled GNSS/INS integration method have changed.

Hence, the system measurement vector of the augmented integration method can be written as:

(30)\begin{equation}{{\bf y}_{\textrm{triple}}}(t) = {[{{\bf y}{{(t)}^T}\textrm{ }{{\bf y}_{\textrm{MM}}}(t)} ]^T}\end{equation}

The new nonlinear measurement model is:

(31)\begin{equation}{{\bf y}_{\textrm{triple}}}(t) = {{\bf h}_{_{\textrm{triple}}}}({\bf x}(t),{{\bf o}_{_{\textrm{triple}}}}(t))\end{equation}

specifically,

(32)\begin{equation}{{\bf h}_{\textrm{triple}}}({\bf x}(t))= [{{\bf h}_\rho }({\bf x}(t))\textrm{ }{{\bf h}_{\dot{\rho }}}({\bf x}(t))\textrm{ }{{\bf h}_{\textrm{MM}}}({\bf x}(t))\textrm{]}\end{equation}

Although the MM position model is linear, it is also conducted using 2n + 1 linear computation as with the pseudorange rate measurement model to simplify the system structure.

The corresponding covariance ${{\bf R}_{\textrm{triple}}}$ of the measurement noise ${{\bf o}_{_{\textrm{triple}}}}(t)$ is up to the GNSS measurement accuracy and DTM accuracy at the same time:

(33)\begin{equation}{{\bf R}_{\textrm{triple}}} = \left[ {\begin{array}{*{20}{c}} {\bf R}& {}\\ {}& {{{\bf R}_{\textrm{MM}}}} \end{array}} \right]\end{equation}

Compared to traditional MM with TC GNSS/INS integration position, the proposed advanced MM measurement-augmented TC GNSS/INS with error-state UKF could, in theory, improve position accuracy in GNSS partly blocked areas. On the one hand, the MM step is carried out in the time update step of the EKF to avoid the poor effect of the degraded GNSS at the current epoch. On the other hand, the matched position as a new measurement to fuse with INS to augment degraded GNSS to generate an accurate position in GNSS partly blocked areas.

Figure 2 displays the architecture of advanced MM measurement-augmented TC GNSS/INS with error-state UKF. The one-step predicted error state output not only for generating sigma points but also for calculating one-step predicted position. Then the DTM is added to help the binary search algorithm to locate the nearest match interval. Based on the nearest match interval, the matched position is generated via a point-line map-matching method. Together, the augmented position measurement and original measurement serve as the input of UKF with different weight.

Figure 2. The architecture of the proposed advanced MM measurement augmented TC GNSS/INS integration with error-state UKF

4.1 Measurement sensing

The GNSS measures the raw pseudorange and pseudorange rate. The IMU measures acceleration and angular velocity to generate the train's position and velocity and is further used to calculate the equivalent pseudorange and pseudorange rate.

4.2 Navigation mode decision-voter

The integration navigation system determines navigation mode using the number of visible satellites. When the train runs in areas where the number of visible satellites is greater than four (defined as GNSS open-sky areas), the system operates in TC GNSS/INS mode with error-state UKF. When the number of visible satellites is fewer than or equal to four and greater than zero (defined as GNSS partly blocked areas), the system switches to advanced MM-augmented TC GNSS/INS mode with error-state UKF. The standalone INS mode is used to provide navigation solutions in GNSS fully blocked areas for a short time. For long tunnels, the MM-augmented INS (INS/MM) mode is used to continue to provide navigation information.

4.3 System filter

Based on this previous decision, the system must estimate the system error state using different system filters.

In TC GNSS/INS mode with error-state UKF, the pseudorange and pseudorange rate difference are input as a measurement. The time updating achieves system error state prediction and corresponding error covariance updating using Equations (11) and (12). Then the UT is used to generate sigma points and associated weights using Equations (13) to (16). The final measurement update is conducted using Equations (19) to (24).

In advanced MM measurement-augmented TC GNSS/INS mode with error-state UKF, the system measurement input includes pseudorange difference, pseudorange rate difference as well as MM position difference. The error-state UKF is still used for the nonlinear estimation.

Based on the different system filters, the proposed robust train localisation method can generate accurate and robust navigation solutions in all environment.

5.1 Navigation performance evaluation of TC GNSS/INS with error-state UKF

To evaluate the navigation performance of the proposed methods, a brief summary of the three methods mentioned in this paper is conducted as follows.

• • Method 1: TC GNSS/INS integration with error-state EKF is abbreviated as EKF GNSS/INS method.

• • Method 2: TC GNSS/INS integration with error-state UKF is abbreviated as UKF GNSS/INS method.

• • Method 3: Proposed advanced MM measurement-augmented GNSS/INS integration with error-state UKF is abbreviated as UKF GNSS/INS/ Advanced MM method.

• • Method 4: Original MM-aided GNSS/INS integration with error-state UKF is abbreviated as UKF GNSS/INS/ Original MM method. Note that the original MM algorithm is detailed in author's previous work (Wei et al., Reference Wei, Sirui, Baigen, Jian and Wei2018).

Figure 7 gives the position difference for method 1 and method 2. The overall position difference of method 2 just has a very slight offset to that of method 1 in all directions. Note that the north position difference for method 2 performs especially smaller compared with method 1. That means the method 2 can generate comparably better accuracy than method 1. In addition to the positioning accuracy, it is obvious that the position solution of method 2 is smoother than that of method 1 in all three directions, especially in areas where visible satellites number fewer than four (1537 ⋅ 3 s to1543 ⋅ 1 s) and satellites frequently change (2069 ⋅ 1 s to 2497 ⋅ 5 s). This is because the posterior mean and covariance of the estimated state with higher accuracy are approximated using a set of weighted sigma points in UKF, which results in a better performance than EKF.

Figure 7. Position difference of method 1 and method 2

Figure 8 plots the velocity difference of method 1and method 2. Similar to the position result, the velocity result of the two methods are also almost coincide to each other except for two jump points.

Figure 8. Velocity difference of method 1 and method 2

To evaluate the navigation performance of method 2 more accurately, Table 1 lists the commonly used statistical results, including position and velocity difference distance root mean Square (DRMS) and root mean square (RMS) of method 1 and method 2. As detailed in Table 1, the position and velocity DRMS and RMS results are almost consistent with the previous analysis, which verifies that UKF can generate equivalent, perhaps even higher-precision, navigation solutions than EKF.

Table 1. Position and velocity comparison of EKF and UKF method

5.2 Robustness performance evaluation of proposed method

The majority of the test environment in this research is GNSS open-sky area, and contains one short GNSS partly blocked area and four short GNSS fully blocked areas. In fact, the practical operation environments are really complex with more and longer GNSS partly blocked scenarios than in the test environment. Since in the areas with severe signal blockage, multipath problems of satellite signals also occur frequently, both the GNSS partly blocked scenario and pseudorange error should be included in the simulation. To further demonstrate the good robustness of the method in GNSS degraded areas, four types of GNSS partly blocked scenarios and one pseudorange error scenario during eight periods were simulated in the test experiments and are listed in Table 2.

Table 2. Simulated scenarios description

Figure 9 shows the position difference comparison between method 1 and method 2 with five types of simulations. For the GNSS partly blocked areas, it is obvious that the position difference of method 1 is larger than that of the method 2 during the seven simulation periods, especially during the only one visible satellite periods (313 s to 333 s), only three visible satellites periods (383 s–483 s and 1800 s to 1900 s) and only four visible satellites periods (2253 s to 2353 s). In the remaining three simulation scenarios, the position difference of the two system performs a light offset to each other. During the pseudorange error periods, the position difference of method 2 is still significantly better that of method 1. This is because the weighting factor in UKF can reduce the impact of pseudorange error on the navigation results.

Figure 9. Position difference of method 1 and method 2 with five types of simulations

The velocity difference for method 1 and method 2 with five types of simulations are also calculated and shown in Figure 10. A similar conclusion can be drawn that the velocity difference of method 2 is smaller and smoother than that of method 1 during the simulation periods.

Figure 10. Velocity difference of method 1 and method 2 with five types of simulations

This position and velocity comparison analysis further confirms that the UKF has better robustness and accuracy than EKF in both GNSS partly blocked areas and pseudorange error areas.

To further observe the robustness of method 2 in GNSS partly blocked areas clearly and concretely, the navigation solutions of method 2 with and without GNSS partly blocked simulations during the seven simulated periods are drawn separately. Due to similar offset tread of method 2 with and without simulations in four simulated scenarios, only one of the typical results in one simulated satellite are analysed in detail, and other results are shown in Appendix A.

Figures 11 and 12 show the position and velocity difference comparison of method 2 with and without one simulated satellite. Note that there is a sight deviation between method 2 with and without one simulated satellite. It is further verified that method 2 can provide robust navigation solutions in GNSS partly blocked environments.

Figure 11. Position difference of method 2 with and without one visible satellite

Figure 12. Velocity difference of method 2 with and without one visible satellite

However, the navigation solution accuracy still has a divergence in different GNSS partly blocked environments and in different navigation directions compared to that in GNSS open-sky environments. Therefore, an advanced map-matching using a binary search algorithm and an additional position update with the DTM is proposed to enhance and optimise navigation performance.

To evaluate the robustness of method 3, the navigation solution comparison of method 2 and method 4 are made. To focus on analysing the robustness only in GNSS partly blocked areas, all the simulated areas are extracted separately and spliced together according to simulation types, as shown in Figure 13. Different simulation scenarios are separated by green dotted lines. The entire simulation time is 570 s.

Figure 13. Comparison of position difference of method 2, method 3 and method 4 in all GNSS partly blocked simulation areas

As shown in Figure 13, the position difference of method 3 in orange, is obviously smaller and smoother than that of method 2 in all three directions in each simulation scenario. It indicates that method 3 can generate better position solutions compared with method 2 in GNSS partly blocked environments. It is because the new MM position measurement generated by advanced map-matching is added to locate with the original GNSS measurements. At the same time, the MM position measurement performs high accuracy and takes up more weight than poor GNSS pseudorange measurement, which is reflected by the corresponding measurement noise covariance matrix.

In addition, the position difference of method 4 is significantly larger than that of method 3 in area 2 with three satellites, but the position difference of both methods is very close to each other in other GNSS partly blocked areas. The north and up position difference of method 4 is larger than that of method 3 in all GNSS partly blocked areas. These results indicate that method 3 can generate better position solutions compared with method 2 in GNSS partly blocked environments. It is because the MM results of advanced MM is not directly output like original MM, but is further integrated with INS as augmented measurement input to weaken the divergence of INS errors.

Hence, the proposed method 3 has superior accuracy and robustness compared with method 2 and method 4 in GNSS partly blocked environments.

Figure 14 shows the comparison of velocity difference of method 2, method 3 and method 4 in all GNSS partly blocked simulation areas. It can be seen that velocity accuracy of method 3 is improved to a certain extent compared with method 2 and method 4. However, compared to the improvement in position accuracy, the improvement in velocity accuracy is not obvious. The MM position measurement is the main reason for the obvious improvement in positioning accuracy. With the aid of the MM position measurement, the GNSS receiver's antenna position changes, which leads to the line-of-sight vector to change, thereby affecting the final velocity estimation.

Figure 14. Comparison of velocity difference of method 2, method 3 and method 4 in all GNSS partly blocked simulation areas

To accurately evaluate positioning accuracy of proposed method 3 in GNSS partly blocked areas, the position and velocity difference DRMS, RMS, mean and standard deviation (STD) of method 2, method 3 and method 4 are summarised in Table 3. The position DRMS and of method 3 is 0 ⋅ 67 m and shows a great improvement of about 71% and 46% compared to that of method 2 in 2 ⋅ 31 m and method 4 in 1 ⋅ 25 m. The up-position RMS also exhibits an improvement of about 29% from 7 ⋅ 69 m in method 2 and about 39% from 8 ⋅ 91 m in method 4 to 5 ⋅ 43 m in method 3.

Table 3. Difference statistics of four different navigation method with simulations

By comparing the position mean, STD and RMS of the two methods, a similar conclusion is obtained. However, the velocity statistical values of the two methods are almost equal.

Furthermore, the statistical results of method 1 (EKF) in GNSS partly blocked areas are also listed in Table 3. The position and velocity DRMS of method 2 have an obvious improvement of about 66% (from 6 ⋅ 75 m to 2 ⋅ 31 m) and 83% (from 0 ⋅ 41 m/s to 0 ⋅ 07 m/s) compared with those of method 1.

The up position and velocity RMS of method 2 also perform an obvious improvement of about 64% (from 21 ⋅ 59 m to 7 ⋅ 69 m) and 81% (from 0 ⋅ 54 m/s to 0 ⋅ 10 m/s). The results further demonstrate that UKF can provide more accurate and robust navigation solutions than EKF when a train runs in GNSS partly blocked areas.