2.1. Principle of residual chi-square detection method

The linear discrete time varying system model can be described as

(1)$$\begin{cases} \boldsymbol{X}_k =\boldsymbol{\varPhi}_{k,k-1} \boldsymbol{X}_{k-1} +\boldsymbol{\varGamma}_{k,k-1} \boldsymbol{W}_{k-1} \\ \boldsymbol{Z}_k =\boldsymbol{H}_k \boldsymbol{X}_k + \boldsymbol{V}_k \\ \end{cases}$$

where $\boldsymbol{X}_{k} $ is the state vector, $\boldsymbol{\varPhi}_{k, k-1} $ is the transition matrix, $\boldsymbol{\varGamma}_{k, k-1} $ is the coefficient matrix, $\boldsymbol{Z}_{k}$ represents the measurement vector, $\boldsymbol{H}_{k}$ is the measurement model matrix. $\boldsymbol{W}_{k}$ is the process noise which is commonly assumed as a zero-mean Gaussian white noise with covariance matrix $\boldsymbol{Q}_{k}$ and $\boldsymbol{V}_{k} $ is the measurement noise which is commonly assumed as a zero-mean Gaussian white noise with covariance matrix $\boldsymbol{R}_{k} $. $\boldsymbol{W}_{k}$ and $\boldsymbol{V}_{k}$ are independent.

One-step prediction of the state in the Kalman filter is

(2)$$\bar{\boldsymbol{X}}_{k/k-1} = \boldsymbol{\varPhi}_{k,k-1} \hat{\boldsymbol{X}}_{k-1}$$

The residual vector is

(3)$$\boldsymbol{v}_k = \boldsymbol{Z}_k - \boldsymbol{H}_k \bar{\boldsymbol{X}}_{k/k-1}$$

The covariance of the residual vector is

(4)$$\boldsymbol{P}_{vk} = \boldsymbol{H}_k \boldsymbol{P}_{k/k-1} \boldsymbol{H}_k^T + \boldsymbol{R}_k$$

When no fault occurs, the residual vector is white noise and its mean is $\boldsymbol{\it 0}$. If the measurement vector contains fault, the statistical characteristics of the residual will change, and its mean is no longer equal to $\boldsymbol{\it 0}$. Define two hypotheses, the null hypothesis H₀ and the alternative hypothesis H₁.

H₀ denotes no fault and can be expressed as

(5)$$\hbox{H}_0 :\boldsymbol{v}_k \sim N (\boldsymbol{\it 0},\boldsymbol{P}_{vk})$$

H₁ denotes the existence of fault and can be given by

(6)$$\hbox{H}_1 : \boldsymbol{v}_k \sim N (\boldsymbol{\mu}, \boldsymbol{P}_{vk})$$

where $\boldsymbol{\mu}$ is the mean of $\boldsymbol{v}_{k}$.

Then, the fault detection function is

(7)$$\Lambda _k = \boldsymbol{v}_k^T \boldsymbol{P}_{vk}^{-1} \boldsymbol{v}_k$$

where Λ_k obeys the χ² distribution and its degree of freedom is the dimension of observations vector $\boldsymbol{Z}_{k}$. Under the null hypothesis, $\Lambda_{k} \sim \chi ^{2} ( 0, \; n) $, under the alternative hypothesis, $\Lambda _{k} \sim \chi^{2} ( \sigma, \; n) $, where σ is the non-centrality parameter and can be obtained as

(8)$$\sigma =\boldsymbol{\mu}^T \boldsymbol{P}_{vk}^{-1} \boldsymbol{\mu}$$

The fault detection criteria are

(9)$$\begin{cases} {\Lambda _k >T_d } & \hbox{the existence of fault}\\ {\Lambda _k \le T_d} & \hbox{no fault occurs} \end{cases}$$

where T _d is the threshold and determines the performance of the fault detection method. According to the Neyman–Pearson criterion, when the false alarm rate is defined as α, the threshold T _d can be worked out through solving the equation $P[ \Lambda _{k} > T_{d} /\hbox{H}_{0}] =\alpha $, and this will minimise the missed detection rate.

2.2. Performance analysis of residual chi-square detection method

Residual chi-square detection is a global detection method for all the filter measurements at time k. To analyse the performance of the residual chi-square detection method, the concept of a minimal detectable bias (MDB) is introduced. The MDB is defined as the model error that can just be detected with the detection method (Yang et al., Reference Yang, Knight, Li and Rizos2013; Teunissen, Reference Teunissen2017).

First, the case in which a fault occurs on only one dimension measurement is considered. If it is assumed that the fault occurs on the ith measurement, then the measurement model can be written as

(10)$$\boldsymbol{Z}_k =\boldsymbol{H}_k \boldsymbol{X}_k +\boldsymbol{V}_k +l_i \nabla$$

where ∇ is the fault in the ith observation, and $l_{i} = ( 0, \; 0, \; \cdots, \; 1, \; \cdots, \; 0) ^{T}$ is the unit vector with the ith element equal to one, when the residual vector is

(11)$$\bar{\boldsymbol{v}}_k = \boldsymbol{Z}_k -\boldsymbol{H}_k \bar{\boldsymbol{X}}_{k/k-1} = \boldsymbol{v}_k +l_i \nabla$$

which has the expectation:

(12)$$E (\bar{\boldsymbol{v}}_k) = \boldsymbol{\mu} =l_i \nabla$$

Then the non-centrality parameter of fault detection function is

(13)$$\sigma _\nabla = \boldsymbol{\mu}^T \boldsymbol{P}_{vk}^{-1} \boldsymbol{\mu} = (l_i \nabla)^T \boldsymbol{P}_{vk}^{-1} l_i \nabla =\nabla ^2 P_{vk \ ii}^{-1}$$

where $P_{vk \ ii}^{-1}$ is the ith diagonal element of covariance matrix $\boldsymbol{P}_{vk}^{-1}$.

Assume the conditional probability density function under the null hypothesis H₀ and the alternative hypothesis H₁ are $p( \Lambda /\hbox{H}_{0}) $ and $p( \Lambda /\hbox{H}_{1}) $. The false alarm rate is defined as α and the missed detection rate is defined as β, the non-centrality parameter σ is the function of α and β, and their relationship can be shown as in Figure 1.

Figure 1. Rates of false alarm and missed detection of residual chi-square method.

From Figure 1, α and β can be expressed as

(14)$$\alpha =\int_{T_d }^\infty {p(\Lambda /\hbox{H}_0)d\Lambda}$$

(15)$$\beta =1-\int_{T_d }^\infty {p({\Lambda /\hbox{H}_1 })d\Lambda } =1-\int_{T_d -\sigma }^\infty {p( {\Lambda /\hbox{H}_0 } )d\Lambda }$$

Combining the above two equations, the non-centrality parameter σ can be given by

(16)$$\sigma =\chi _\alpha ^2 (n)-\chi _{1-\beta }^2 (n)$$

The MDB for the ith observation can then be obtained according to Equations (13) and (16).

(17)$$| \nabla|_{\min } =\sqrt {\sigma /P_{vk \ ii}^{-1}}$$

If a fault occurs on multiple dimensions of measurements at the same time, assume that the fault vector is $\nabla = ( {\nabla _{1}, \; \nabla _{2}, \; \cdots, \; \nabla _{i}, \; \cdots, \; \nabla _{n} }) ^{T}$, where ∇_i is the fault value in the ith observation; if there is a fault, $\nabla_{i} \ne\ 0$, otherwise ∇_i = 0. The corresponding coefficient matrix is $\boldsymbol{L}$. Under this case, the non-centrality parameter σ can be written as

(18)$$\bar {\sigma }_\nabla = \boldsymbol{\mu}^T \boldsymbol{P}_{vk}^{-1} \boldsymbol{\mu} = ( {\boldsymbol{L}\nabla } )^T \boldsymbol{P}_{vk}^{-1} \boldsymbol{L}\boldsymbol{\nabla}$$

Recalling Equation (16), a fault can be detected easily only when the fault vector meets the following equation:

(19)$$\sigma \le ({\boldsymbol{L}\boldsymbol{\nabla} })^T \boldsymbol{P}_{vk}^{-1} \boldsymbol{L}\boldsymbol{\nabla}$$

The larger the fault, the more easily it can be detected. When the fault is smaller than the MDB, it is difficult to detect the fault accurately and the probability of missed detection will increase. Faults do not always have a large size of value in practical applications, so there is a need to make an improvement in the residual chi-square detection method.

Moreover, the residual chi-square detection method is a system-level evaluation. After detection of the fault, the faulty subsystem is isolated and the state estimation is obtained by using the one-step prediction of the Kalman filter in general. On the one hand, isolating the faulty subsystem is a waste of useful measurements and, as a result, the estimation accuracy will decrease. On the other hand, the error of the one-step prediction of the Kalman filter will increase with the recursive time, and this in turn will have an influence on the value of the residual and the fault detection sensitivity will degrade, which would cause more of the missed detection phenomenon. To sum up, the performance of the residual chi-square detection method is affected by the size and duration time of the fault. To improve the fault detection performance, three improvements are introduced and the new fault detection scheme is designed in the next section.

3.1. Fault detection and identification based on local test

The residual chi-square detection method is a system-level method, it cannot identify the faulty measurement. But identifying and isolating the faulty measurement in time is of great significance to improve the performance of fault detection and the precision of the integrated navigation system. To identify the faulty measurement and retain more useful normal observations, the local test method is introduced.

Assume the fault vector at time k is $\boldsymbol{\nabla}_{k} = ( \nabla _{k}^{1}, \; \nabla _{k}^{2}, \; \cdots, \; \nabla _{k}^{i}, \; \cdots, \; \nabla _{k}^{n}) $, the residual vector is $\bar{\boldsymbol{v}}_{k}$, then the fault detection function can be given as

(20)$$\lambda _k^i =\frac{\bar {\boldsymbol{v}}_k^i }{\sqrt {\boldsymbol{P}_{vk}^{ii} } }$$

where $\boldsymbol{P}_{vk}^{ii}$ is the ith diagonal element of covariance matrix $\boldsymbol{P}_{vk}$, and $\bar {\boldsymbol{v}}_{k}^{i}$ is the ith element of $\bar{\boldsymbol{v}}_{k}$. $\lambda _{k}^{i} $ is the standardised residual.

If there is no fault on the ith dimension of the measurement vector, then the null hypothesis is $\hbox{H}_{0} \colon \lambda _{k}^{i} \sim N( {0, \; 1}) $, otherwise the alternative hypothesis is $\mbox{H}_{1} \colon \lambda _{k}^{i} \sim N( {\delta _{i}, \; 1}) $, where δ_i is the non-centrality parameter, and it is the function of the false alarm rate α₀ and missed detection rate β₀. Their relationship can be shown in Figure 2. For the sake of writing the later equations conveniently, we replace $\lambda _{k}^{i} $ with its absolute value, and $\lambda _{k}^{i} $ in the later part is equal to $\vert {\bar {\boldsymbol{v}}_{k}^{i} } \vert /\sqrt {\boldsymbol{P}_{vk}^{ii}}$.

Figure 2. Rates of false alarm and missed detection of local test method.

From Figure 2, the non-centrality parameter can be expressed as

(21)$$\delta =N_{1-\alpha _0 /2} (0,1)-N_{\beta _0 } (0,1)$$

The fault detection criteria are

(22)$$\begin{cases} {\lambda _k^i >N_{\alpha _0 /2} ( {0,1})} & \hbox{fault occurs}\\ {\lambda _k^i \le N_{\alpha _0 /2} ( {0,1})} & \hbox{no fault occurs}\\ \end{cases}$$

Ideally, the probability of false alarm in the ith pseudorange α₀ is set so that the probability of any one of the outlier tests failing, when the null hypothesis is true, is equal to the global test probability of false alarm α (Yang et al., Reference Yang, Knight, Li and Rizos2013). Usually, the relationship of α₀ and α can be expressed as in Kelly (Reference Kelly1998):

(23)$$\alpha _0 =1-\sqrt[n]{1-\alpha }$$

In the same way, the missed detection rate of local test β₀ can be approximated as β (Kelly, Reference Kelly1998). The MDB for the ith observation can be obtained as

(24)$$| {\nabla _i } |_{\min } =\delta \sqrt {\boldsymbol{P}_{vk}^{ii} }$$

Every dimensional measurement should be tested to identify the fault because a fault may occur on each measurement. After identification of the fault, the faulty measurement is isolated, while the remaining normal measurements are used to conduct the filter measurement update process.

3.2. Fault detection based on robust estimation

In practical application, the value of a fault is not always large; small faults and soft faults often occur, and these types of fault are hard to detect. The traditional local detection method gives the alarm when the detection function is larger than the threshold, otherwise it indicates no fault. According to the analysis in Subsection 3.1, when the measurement contains a small size of fault, the detection method can easily give the missed detection result. That is to say, when the test function is smaller than the threshold, it may include two cases: one is that no fault occurs on the measurement in reality, the other is that the fault is too small to detect, and this may result in decrease of filter accuracy and fault detection sensitivity. Hence, in order to improve the performance of the local test method for the small size of fault, the robust estimation scheme is introduced. The weight function of the measurements can be defined as follows:

(25)$$p_i = \begin{cases} 1 & {\lambda _k^i \le T_{d1} } \\ {\dfrac{T_{d1} }{\lambda _k^i } \left( {\dfrac{T_{d2} -\lambda _k^i }{T_{d2} -T_{d1} }}\right)^2} & {T_{d1} <\lambda _k^i \le T_{d2} } \\ 0 & {T_{d2} <\lambda _k^i }\end{cases}$$

where T _d1 and T _d2 are constant and selected as $T_{d1} =1\cdot 0-1\cdot 5$ and $T_{d2} =2\cdot5-8$ in general (Yang et al., Reference Yang, He and Xu2001). It can be seen from Equation (25) that when the fault detection function $\lambda _{k}^{i} $ is smaller than T _d1, it can be considered that no fault occurs on the ith dimensional measurement, the weight of this dimensional measurement will be equal to 1. When $\lambda _{k}^{i} $ is larger than T _d2, it indicates that there is a fault on the ith dimensional measurement and the weight of this measurement will be selected as 0 to isolate its influence on the filter precision. When $\lambda _{k}^{i} $ is between T _d1 and T _d2, it means that the ith dimensional measurement may be polluted by fault error and the quality of this measurement is lower than the normal measurement, so a weight factor smaller than 1 will be used to reduce the influence of this measurement on the filter precision. Therefore, T _d1 and T _d2 can be regarded as two thresholds for faults with different values.

According to the adjustment of the weight function, each dimension of the measurements will be given different weight based on the quality of the measurement. And this will affect their weight in the filtering estimation. The better the quality of the measurement, the greater the weight.

Define $\bar {\boldsymbol{P}}_{k} $ as the equivalent weight matrix, and

(26)$$\bar {P}_k^i =R_k^i /p_i$$

where p _i is the nonzero element, $\bar {P}_{k}^{i} $ and $R_{k}^{i} $ are the ith diagonal elements of $\bar {\boldsymbol{P}}_{k} $ and $\boldsymbol{R}_{k} $, respectively. Then the new filter gain matrix can be written as

(27)$$\hat {\boldsymbol{K}}_k =\boldsymbol{P}_{k/k-1} \boldsymbol{H}_k^T ( {\boldsymbol{H}_k \boldsymbol{P}_{k/k-1} \boldsymbol{H}_k^T +\bar {\boldsymbol{P}}_k })^{-1}$$

When a small fault occurs but the fault detection function is lower than T _d2, the weight of the faulty measurement will be adjusted to be less than 1. This means that the weight of this measurement in the Kalman filter measurement update will decrease and the filter accuracy will be improved, in return this would improve the fault detection performance the next time.

According to Equation (25), it can be concluded that T _d1 and T _d2 will have a great influence on the performance of the filter, so suitable thresholds need to be selected. It can be seen from Equation (25) that when the fault detection function $\lambda _{k}^{i} $ is larger than T _d2, the ith dimensional measurement will be isolated, so T _d2 is taken as the threshold, which is equal to the threshold $N_{\alpha _{0} /2} ( {0, \; 1}) $ of the local test.

In the next part, the effect of T _d1 on the performance of filter will be analysed. Define

(28)$$f ( {T_{d1} })=\frac{T_{d1} }{\lambda _k^i } \left( {\frac{T_{d2} -\lambda_k^i }{T_{d2} -T_{d1} }} \right)^2$$

For a fixed $\lambda _{k}^{i} $ and T _d2, f(T _d1 ) can be seen as the function of T _d1. Take the derivative of f(T _d1 ), with respect to T _d1, the derived function can be written as

(29)$${f}'( {T_{d1} } )=\frac{( {T_{d2} -\lambda _k^i } )^2( {T_{d2} +T_{d1} } )}{\lambda _k^i ( {T_{d2} -T_{d1} } )^3}$$

From the above equation, it can be concluded that f′(T _d1 ) > 0 when $T_{d1} \in ( {0, \; T_{d2} } ) $, which means f(T _d1 ) is a monotone increasing function when T _d1 is in the interval (0, T _d2 ). In other words, the weight of the measurement will increase with the value of T _d1. Assume that $\lambda _{k}^{i} $ is the standardised residual of the ith measurement, and $T_{d1} < \lambda _{k}^{i} < T_{d2} $. If a fault occurs on this measurement, the larger the value of T _d1, the greater the weight of the faulty measurement, and the worse the precision of navigation. In contrast, if no fault occurs, the larger the value of T _d1, the greater the weight of the normal measurement, and the better the precision of navigation. Therefore, the value of T _d1 will have a great effect on the weight of both normal and faulty measurements. The optimal value is difficult to choose, and usually the choice of the T _d1 is based on the experience or the requirement of the practical application.

As the probability of the fault is much smaller than that of normal operation, the performance of filtering with different values of T _d1 is investigated when no fault occurs. The results show that the positioning error will increase by nearly 30% and 100% when T _d1 = 1.5 and T _d1 = 1, respectively. And the false alarm rate corresponding to T _d1 = 1.5 is more than 100 times the false alarm rate α₀ of the local test. When the fault detection function $\lambda _{k}^{i} $ is smaller than 1.5, the ith dimensional measurement can be regarded as normal measurement. Based on the above analysis, to improve the fault detection performance for small size of fault and ensure the weight of normal measurements, T _d1 is selected to be equal to 1.5.

3.3. Improved fault detection method aided by sliding window test

In Subsection 3.2, when the fault detection function $\lambda _{k}^{i} $ is between T _d1 and T _d2, a weight factor is added to reduce the weight of the ith observation. In fact, a small size of fault may have occurred but the local test based on the standardised residual fails to detect it. Restricted by the detection ability of the local test, the performance of the fault detection for small faults is still not ideal and remains to be further improved. Therefore, to improve the detection performance for small faults, a sliding window test based on the past measurements is added.

In the sliding window test, the measurements used are not limited to the current epoch. Assume that the window width is m, the new residual vector is defined as

(30)$$\hat {\boldsymbol{v}}_k = \left( {\hat {\boldsymbol{P}}_{vk}^{-1} }\right)^{-1} \mathop{\boldsymbol{\sum}}\nolimits_{j=k-m+1}^k {\boldsymbol{P}_{vj}^{-1} \boldsymbol{v}_j }$$

where $\hat {\boldsymbol{P}}_{vk} $ and $\boldsymbol{P}_{vj} $ are the covariance matrix of $\hat{\boldsymbol{v}}_{k} $ and $\boldsymbol{v}_{j} $, respectively. $\hat {\boldsymbol{P}}_{vk}$ can be expressed as

(31)$$\hat {\boldsymbol{P}}_{vk}^{-1} = \mathop{\boldsymbol{\sum}}\limits_{j=k-m+1}^k {P_{vj}^{-1} }$$

From Equations (30) and (31), it can be seen that the new residual vector of the sliding window test is a weighted average of the Kalman filter residual over the past measurements, and it can be obtained that

(32)$$\hat{\boldsymbol{P}}_{vk}^{ii} <\boldsymbol{P}_{vk}^{ii}$$

where $\hat{\boldsymbol{P}}_{vk}^{ii}$ and $\boldsymbol{P}_{vk}^{ii}$ are the ith diagonal element of $\hat{\boldsymbol{P}}_{vk}^{ii}$ and $\boldsymbol{P}_{vk}$, respectively. The mathematical derivation of Equation (32) is as follows.

Proof. $\hat{\boldsymbol{P}}_{vk} $ and $\boldsymbol{P}_{vj}$ are the covariance matrix of $\hat{\boldsymbol{v}} _{k}$ and v _j, so $\hat{\boldsymbol{P}}_{vk} $ and $\boldsymbol{P}_{vj}$ are symmetric positive definite matrices. When the window width m is 2, Equation (31) can be written as

(33)$$\hat{\boldsymbol{P}}_{vk}^{-1} =\boldsymbol{P}_{vk-1}^{-1} +\boldsymbol{P}_{vk}^{-1}$$

For the sake of convenience in writing, Equation (33) can be written as

(34)$$\boldsymbol{P}_g^{-1} =\boldsymbol{P}_1^{-1} +\boldsymbol{P}_2^{-1}$$

where $\boldsymbol{P}_{g} =\hat{\boldsymbol{P}}_{vk} $, $\boldsymbol{P}_{1} =\boldsymbol{P}_{vk-1} $, $\boldsymbol{P}_{2} =\boldsymbol{P}_{vk} $. Based on the matrix theory, Equation (34) can be written as

(35)$$\begin{align} \boldsymbol{P}_g &=(\boldsymbol{P}_1^{-1} +\boldsymbol{P}_2^{-1} )^{-1} \\ &=\boldsymbol{P}_1 ( {\boldsymbol{P}_1 +\boldsymbol{P}_2 } ) ^{-1}\boldsymbol{P}_2 \\ & =\boldsymbol{P}_1 -\boldsymbol{P}_1 ( {\boldsymbol{P}_1 +\boldsymbol{P}_2 } ) ^{-1}\boldsymbol{P}_1 \\ &=\boldsymbol{P}_1 -\boldsymbol{P}_1^T ( {\boldsymbol{P}_1 +\boldsymbol{P}_2 } ) ^{-1}\boldsymbol{P}_1 \\ & =\boldsymbol{P}_1 -\hat {\boldsymbol{P}}_1 \end{align}$$

Because $\boldsymbol{P}_{1} $ and $\boldsymbol{P}_{2} $ are symmetric positive definite matrices, $\boldsymbol{P}_{1} +\boldsymbol{P}_{2} $ is also a symmetric positive definite matrix; therefore, the following equation can be obtained:

(36)$$\label{eq36} \hat{\boldsymbol{P}}_1^{ii} =\boldsymbol{P}_{1i}^T ( {\boldsymbol{P}_1 +\boldsymbol{P}_2 } ) ^{-1}\boldsymbol{P}_{1i} >0$$

where $\hat{\boldsymbol{P}}_{1}^{ii} $ is the ith diagonal element of $\hat{\boldsymbol{P}}_{1} $, and $\boldsymbol{P}_{1i} $ is the ith column of $\boldsymbol{P}_{1} $. Then, it can be obtained that the following equation is established:

(37)$$\boldsymbol{P}_g^{ii} <\boldsymbol{P}_1^{ii}$$

where $\boldsymbol{P}_{g}^{ii} $ and $\boldsymbol{P}_{1}^{ii} $ are the ith diagonal element of $\boldsymbol{P}_{g} $ and $\boldsymbol{P}_{1} $, respectively.

In the same way, the following equation can be obtained:

(38)$$\boldsymbol{P}_g^{ii} < \boldsymbol{P}_2^{ii}$$

where $\boldsymbol{P}_{2}^{ii} $ is the ith diagonal element of $\boldsymbol{P}_{2} $. So it can be concluded that $\boldsymbol{P}_{g}^{ii} < \boldsymbol{P}_{l}^{ii}, \; \exists l\in ( 1, \; 2) $. The same result can be easily obtained when m > 2 by using mathematical induction. Therefore, Equation (32) is proved.

The test statistic is given by

(39)$$\hat {\lambda }_k^i =\frac{\vert {\hat{\boldsymbol{v}}_k^i } \vert }{\sqrt {\hat {\boldsymbol{P}}_{vk}^{ii} } }$$

The fault detection criteria are the same as in Equation (22).

According to the analysis in Subsection 3.1, the MDB of the ith observation in this method can be obtained as

(40)$$\vert {\hat {\nabla }_i } \vert _{\min } =\delta \sqrt{\hat {\boldsymbol{P}}_{vk}^{ii}}$$

Combining Equations (24), (32) and (40), it can be obtained that

(41)$$\vert {\hat {\nabla }_i } \vert _{\min } <\vert {\nabla_i } \vert _{\min }$$

Equation (41) indicates that the MDB of the sliding window test is smaller than that of the local test based on current measurement. In other words, the sliding window test has better sensitivity for the small size of fault than the local test based on current measurement. But when no fault occurs, due to the higher sensitivity of fault detection, the fault detection function of the sliding window test may be greater than that of the local test based on current measurement, which may result in more normal measurements being judged as polluted or faulty measurements. If using the sliding window test in the whole detection, the weight of more normal measurements may reduce, which will result in a decrease of the navigation precision when no fault occurs. Besides, it also involves a heavy computation burden due to the storage of historical data and the calculation of covariance matrix inversion. On this basis, the sliding window test is added to assist the fault detection method in Subsection 3.2 as follows.

For i = 1, 2,···, n, the fault detection function of the ith dimensional measurement at time k − 1 is $\lambda _{k-1}^{i} $, and a fault may have occurred on the ith observation when

(42)$${T_{d1} < \lambda _{k-1}^i \le T_{d2} } \quad {\exists i\in ({1,2,\cdot \cdot \cdot, n})}$$

To improve the fault detection performance, the sliding window test is adopted to conduct the fault detection process at time k. If the above equation is not established, then the local test based on the current measurement is used in the fault detection process of epoch k. By using this strategy, the local test is adaptively adjusted between the sliding window test and the traditional local test based on current measurements.

The fault detection scheme of this paper is shown in Figure 3. According to Figure 3, the steps of the improved fault detection method are as follows.

1. (1) Calculation of the residual vector. Choose the fault detection method according to the fault detection function of the last epoch. If Equation (42) is satisfied, the sliding window test is adopted and the residual is calculated by Equation (30), otherwise, the local test based on the current measurement is used and the residual vector is calculated by Equation (11).

2. (2) Standardisation of the residual. Calculate the standardised residual of each dimensional measurement. If the residual is calculated by Equation (30), the standardised residual can be obtained by $\hat {\boldsymbol{v}}_{k}^{i} /\sqrt {\hat {\boldsymbol{P}}_{vk}^{ii} } $. On the contrary, if the residual is calculated by Equation (11), the standardised residual will be calculated by Equation (20).

3. (3) Local fault detection and identification. Compare the standardised residual $\lambda _{k}^{i} $ with the two thresholds T _d1 and T _d2. If $\lambda _{k}^{i} \le T_{d1} $, it indicates that no fault occurs on the ith dimension measurement, the weight factor is 1. If $T_{d2} < \lambda _{k}^{i} $, there is fault on the ith dimension measurement, the weight factor is 0 and the measurement is isolated. If $T_{d1} < \lambda _{k}^{i} \le T_{d2} $, calculate the weight factor according to Equation (25).

4. (4) Filter measurement update. After all the local tests are complete, calculate the equivalent weight matrix, and then conduct the filter measurement update.

5. (5) Fault detection process at time k is complete, return to step 1 at time k + 1.

Figure 3. The fault detection scheme.

4.1. Simulation conditions

The gyro constant drift is 0.1 h with its random drift 0.1 h, the constant bias and random bias of accelerometer is 50 μg. The initial geographical position is 108^° east longitude and 34^° north latitude, the initial velocity is zero, and the initial azimuth is 90^°. The flight trajectory includes acceleration, deceleration, climbing motion, diving motion and turning motion. The number of visible satellites is eight (S1 to S8) and the pseudorange measurement error standard deviation is 10 m. The output frequency of the inertial measurement unit (IMU) and GNSS receiver is 100 Hz and 1 Hz, respectively. The whole simulation time is 1,600 s, and the filter cycle is set to be 1 s.

The state vector is $ \boldsymbol{X}=[\phi _{E}, \; \phi _{N}, \; \phi _{U}, \; \delta v_{E}, \; \delta v_{N}, \; \delta v_{U}, \; \delta L, \; \delta \lambda, \; \delta h, \; \varepsilon _{x}, \; \varepsilon _{y}, \; \varepsilon _{z}, \; \nabla _{x}, \; \nabla _{y}, \; \nabla _{z}, \; \delta t_{u}, \; \delta t_{ru} ]^{T}$, where ϕ _E, ϕ _N, ϕ _U are the misalignment angles, δ v _E, δ v _N and δ v _U are the east, north and upward velocity errors, respectively. δ L, $\delta \lambda $ and δ h denote the latitude, longitude and height errors. $\varepsilon _{x} $, $\varepsilon _{y} $, $\varepsilon _{z} $ and ∇_x, ∇_y, ∇_z represent the gyro biases and accelerometer biases. δ t _u, δ t _ru are the range bias and range drift related to the receiver clock. The observations of the filter are the pseudorange differences between INS and GNSS.

To verify the performance of the proposed method (M4), the other three methods are introduced in the simulations for comparison. The first one is the residual chi-square detection method (M1), the second is the local fault detection method in Subsection 3.1 (M2), the third is the fault detection method proposed in Subsection 3.2 (M3).

The false alarm rate of M1 is 0.01 and its missed detection rate is 0.2. So it can be obtained by calculation that the thresholds of M1 and M2 are 20.09 and 3.226, respectively. For M3 and M4, the thresholds T _d1 and T _d2 are selected as 1.5 and 3.226, according to the analysis in Subsection 3.2.

4.2. Simulation results

First the MDB of the residual chi-square detection method and the local fault detection method is analysed. Assume that only one dimensional measurement has fault, then the MDB of M1 and M2 for each measurement is shown in Figure 4.

Figure 4. MDB of M1 and M2 for each measurement.

Under the case that only one measurement has fault, from Figure 4 it can be seen that the MDB of each measurement of M1 is about 39.5, while that of M2 is about 40.8, which is a little larger than that of M1. As the residual chi-square detection method is a system-level method, the fault detection result is the comprehensive effect of the faults on all the measurements. In other words, if small faults occur on multiple dimensional measurements at the same time, it may result in the situation that the fault can be detected by M1 but M2 fails to detect the fault.

To choose an appropriate window width, the MDB of the sliding window test with different window widths is shown in Figure 5. All the MDB are calculated from time k = 100. From Figure 5, it can be seen that the larger the window width, the smaller the MDB. In other words, the fault detection sensitivity for a small fault will increase with the enlargement of the window width. The MDB for three different window widths are about 18 m, 13 m and 9 m, respectively. Although the sensitivity of fault detection will increase with the enlargement of the window width, it will also bring greater calculation burden. Besides, the enlargement of the window width will further reduce the weight of the normal measurement when no fault occurs and result in the decrease of navigation precision. Therefore, to give the whole navigation process high navigation precision and to impose less computation burden, the window width is chosen as 5 m.

Figure 5. MDB of the sliding window test with different window widths.

To verify the fault detection performance of the proposed method, several simulation tests with different fault models are performed.

Test 1: Assume that a large abrupt fault occurs in the fourth satellite pseudorange observation in the period of 600–900 s, and the value of the fault is 100 m. The detection results of the four methods are shown in Figure 6. Here, for the convenience of comparison, only the fault detection function of the fourth satellite pseudorange is given in M2, M3 and M4.

Figure 6. Fault detection function of four methods: Test 1, large abrupt fault. (a) M1 and M2; (b) M3 and M4.

From Figure 6, it can be concluded that the values of the fault detection function of four methods obviously increase when the abrupt fault occurs at 600 s. And during the interval of 600–900 s, the values of the fault detection function are larger than the thresholds, which indicates that the large abrupt fault can be detected by all the four methods. In Figure 6(a) the fault detection function of M1 has a trend of decline. It can be imagined that if the fault interval lasts for a long time, the fault detection function may decrease to lower than the threshold, which would lead to missed detection. The fault detection functions of M2 and M3 fluctuate near a stable value, and the fault detection performance is not affected by the duration of the fault. From Figure 6(b), it can be seen that the fault detection function values often jump due to the frequent switching between the local test based on the current measurement and the sliding window test based on past measurements.

Test 2: Assume that a small abrupt fault occurs in the fourth satellite pseudorange in the period of 600–900 s, and the value of the fault is 30 m. The detection results of the four methods are shown in Figure 7. Only the fault detection function of the fourth satellite pseudorange is given in M2, M3 and M4.

Figure 7. Fault detection function of four methods: Test 2, small abrupt fault. (a) M1 and M2; (b) M3 and M4.

From Figure 7, it can be seen that a large number of the fault detection function values for the four methods are lower than the threshold during the fault stage when a small fault occurs. That means that a lot of missed detection phenomena appear. Compared with M1 and M2, M3 adopts a robust filter and it can reduce the weight of the faulty measurement when it fails to detect and isolate it, so it has better fault detection performance and less missed detection phenomena. But from Figure 7(b), it can be seen that there are still a large number of detection function values which are between the two thresholds or under threshold T _d1 in M3. M4 adopts the sliding window test which has better detection performance than the local test which only uses the current measurement, therefore, M4 has less missed detection phenomena compared with M3.

Test 3: Assume that multiple measurements have faults at the same time. Here, the 30 m, 45 m and 60 m constant errors are added to the fourth, sixth and eighth satellite measurements. Again, the fault stage is from 600  s to 900 s. The detection results of the four methods are shown in Figure 8. Here, the fault detection functions of S4, S6 and S8 are given for M2, M3 and M4.

Figure 8. Fault detection function of four methods: Test 3, multiple faults at the same time. (a) M1; (b) M2, M3 and M4 for the fourth, sixth and eighth satellites.

It can be seen from Figure 8(a) that when different size faults occur on several measurements at the same time, although the value of the fault on each measurement is not very large, as the combined effect of faults on all the measurements is large enough, M1 can still identify the fault accurately. Figure 8(b) shows the fault detection function of the other three methods. It can be concluded that all three methods have a lot of missed detection phenomena for small faults, and the larger the fault, the less the missed detection phenomena. But for M2, the detection function value of S6 is larger than threshold for a long time after the fault disappears at 900 s, which results in a large number of false alarms. Compared with M2, M3 and M4 has better performance and less missed detection phenomena, and M4 has the best detection performance.

Test 4: The above three tests all aimed to verify the fault detection performance of the proposed method regarding abrupt faults. In practical applications, as well as abrupt faults, gradual faults are of common occurrence and difficult to detect. So, in order to prove the validity and superiority of the proposed method in detecting gradual faults, assume that a gradual fault whose slope is 0·5 m/s occurs on the fourth satellite measurement in the period of 600 ~ 800 s. The detection results of the four methods are shown in Figure 9.

Figure 9. Fault detection function of four methods: Test 4, gradual fault. (a) M1 and M2; (b) M3 and M4.

It can be seen from Figure 9 that it takes a long time for all the four methods to detect the gradual fault. The fault detection time delay of each method is about 110 s, 100 s, 70 s and 60 s, respectively. This means that the proposed method can detect the gradual fault the fastest and has less missed detection phenomena. The reason is that the influence of undetected faults on filter estimation is reduced in M4, whose fault detection sensitivity is better than those of M1, M2 and M3.

The false alarm rate and missed detection rate are usually used to measure the performance of the fault detection method, and the navigation precision is usually measured by the root mean square error (RMSE). To further illustrate the performance of the proposed fault detection method, a 100-run Monte Carlo simulation of the four methods in each test was conducted to obtain statistical results. The comparison of fault detection results and the root mean square of horizontal position errors of the four methods under four tests are shown in Table 1 and Figure 10, respectively.

Figure 10. RMSE of horizontal position errors of four methods in the four tests. (a) Test 1; (b) Test 2; (c) Test 3; (d) Test 4.

Table 1. Comparison of fault detection results of four methods under four tests.

In test 3, faults occur on three satellites' measurements, so the fault detection results of S4, S6 and S8 are listed in Table 1 for M2, M3 and M4, while the global fault detection results are given for M1 as a result of being a system-level method. From Table 1, it can be seen that all the four methods have good fault detection performance when the value of the fault is large enough. When a small size of fault occurs, the performance of the traditional residual chi-square detection method and the local test is poor, and the missed detection rate is higher than 80%. With the aid of robust estimation, the missed detection rate will be reduced by 16%, but there is still a large number of missed detection phenomena. The proposed method can effectively improve the fault detection performance, and the missed detection rate can be reduced by 49% compared with the conventional method. For the false alarm rate, the given theoretical value of M1 is 0.01, while that of the other three methods can be calculated as 0.0013 according to Equation (23). It can be seen from Table 1 that the false alarm rate of the proposed method is higher than those of the traditional methods due to the higher sensibility of fault detection, but the value is tolerated and the increase of false alarm rate has little impact on the continuity and reliability of the integration. Therefore, the proposed method has better performance of fault detection.

It can be seen from Figure 10 that the proposed method M4 has high precision under all four tests. When a large size of fault occurs, although the faulty measurement is isolated by M1, the position error would accumulate till the fault disappears because the error of one-step prediction increases with time. For M2, M3 and M4, the faulty measurements are isolated and other normal measurements are used for the state estimation, therefore, they have higher precision. When a small fault or gradual fault occurs, as there is a high rate of missed detection for M1, M2 and M3, the position error is still rather large, while M4 has better fault detection and tolerance performance, so it has higher precision.

The mean RMSEs of the four methods during the fault period are shown in Table 2. It can be seen that when a small fault or a gradual fault occurs, the proposed method has higher navigation accuracy compared with the other three methods, and the positioning accuracy is obviously superior to the traditional methods M1 and M2. When the value of the fault is large enough, as in test 1, for M3, as the result of adopting the robust estimation, the weight of the normal measurement will be reduced, therefore, the navigation accuracy is a little lower than the conventional method M2. For the proposed method M4, because of combining the robust estimation and sliding window test, the weight of normal measurement will be further reduced, so the position precision of M4 is a little lower than that of M3. Besides, when multiple faults occur in test 3, although the value of the fault on each measurement is not very large, M1 still has good performance at fault detection due to being a global detection method. However, M2 is a local test method, and the faulty measurements are difficult to detect because the value of fault is small, which results in missed detection and the decrease of the filtering precision, so the RMSE of M2 is larger than that of M1. Compared with M2, M3 and M4 have better performance of fault detection and tolerance for small size of fault because of adopting the robust estimation method. Therefore, they have higher filtering precision than M2. For M3, as a result of the large number of missed detection phenomena, the performance of fault detection and tolerance for small size of fault is still a little poor, so it has a larger RMSE than M1.

Table 2. Comparison of mean RMSE of horizontal position errors during fault period.

In practical application, the uncertainty of the fault model will have a great influence on the performance of the fault detection method. From the above simulation results, it can be seen that the proposed method has better adaptability for the fault model and this would be of great importance to improve the reliability and robustness of the navigation system.