A Novel Method of Fault Detection and Identification in a Tightly Coupled, INS/GNSS-Integrated System

Fault detection and identification are vital for guaranteeing the precision and reliability of tightly coupled inertial navigation system (INS)/global navigation satellite system (GNSS)-integrated navigation systems. A variance shift outlier model (VSOM) was employed to detect faults in the raw pseudo-range data in this paper. The measurements were partially excluded or included in the estimation process depending on the size of the associated shift in the variance. As an objective measure, likelihood ratio and score test statistics were used to determine whether the measurements inflated variance and were deemed to be faulty. The VSOM is appealing because the down-weighting of faulty measurements with the proper weighting factors in the analysis automatically becomes part of the estimation procedure instead of deletion. A parametric bootstrap procedure for significance assessment and multiple testing to identify faults in the VSOM is proposed. The results show that VSOM was validated through field tests, and it works well when single or multiple faults exist in GNSS measurements.


Introduction
The inertial navigation system (INS) and global navigation satellite system (GNSS) compose the tightly coupled integrated navigation system directly. The fusion of the raw navigation information (pseudo-range or carrier phase measurements) of the GNSS and the inertial measurements of IMUs (inertial measurement units) is implemented in some nonlinear filters [1,2], such as the cubature Kalman filter [3] and the unscented Kalman filter [4]. An INS is a self-contained dead-reckoning system that does not rely on external information and is immune to interference. It can be assumed to be perfectly reliable. GNSS measurements are more vulnerable and are more easily jammed or interfered with. GNSS measurements may be disturbed by faults [5]. The integrated system provides superior performance when compared to either a stand-alone INS or GNSS due to their complementary characteristics [6]. Fault detection and identification play important roles in the practical applications. If a fault is not detected and identified instantly, the navigation solution will be contaminated by the fault in the measurements, and the precision and reliability will degrade. Therefore, the need for fault detection and identification in integrated navigation systems is paramount [7].
Two main types of fault detection and identification are usually employed in tightly coupled integrated navigation systems. The most common methods are fault detection and isolation (FDI), fault detection and exclusion (FDE) and fault detection and recovery (FDR). When a fault of the raw sensor signal in the integrated navigation system occurs, the FDE can provide an alarm and enable the navigation system to exclude the faulty measurements. A fault detection algorithm based on hypothesis testing in parity space was investigated by Sturza [8]. Integrity and quality control can be implemented through recursive filtering and residual chi-square tests [9]. Detection of the state chi-square relying on the residual of the current epoch can detect the abrupt changing faults, but the method does not deal with the detection of gradual changing faults well. Based on hypothesis testing, autonomous integrity monitoring by an extrapolation method was investigated for detecting gradual faults, and the measurements used in this method derive not only from current epoch but also from the previous epochs [10]. Extended receiver autonomous integrity monitoring introduced the error model of the nonlinear filter into the monitoring process [11]. Two independent detectors exist for GNSS faults and filter faults: an exclusion function can be utilized for the identification and removal of the faulty measurements, and elimination of the filter fault effect is carried out through filter recovery [12]. The former method involves an adaptive filter or a robust filter.
An adaptive filter guarantees the precision and reliability of the integrated navigation system through adaptive adjustments of the noise covariance matrix and reduction of the weights of faulty measurements, and avoids hypothesis testing. There are various kinds of adaptive estimation methods, such as the generalized maximum likelihood estimator [13], the soft-threshold optimal estimator [14] and median least squares [15]. A robust filter can reduce the weights of the fault measurements in the estimation, and its performance mainly rests on the selection of the weight matrix. Crespillo proposed robust M-estimators [16]. Compared with classical extended Kalman filters, robust M-estimators offer increased resilience at the estimator level and limit the different faulty effects in the final estimation. Appropriate thresholds are utilized for calculating the weight factor for each measurement, and adjusting the gain matrix adaptively to reduce the influence of the undetected faulty measurement. The sliding window test contributes to the improvement of the fault detection performance [17]. The filters do not need to delete the faulty measurements.
Artificial intelligence has been applied to enhance the performance of fault detection, and benefits from rapid development. A data-driven adaptive neuron fuzzy inference system for predefined faults was used for the detection of navigation sensor faults in unmanned aerial vehicles [18]. Gaussian process regression was first utilized to calculate the innovation of a Kalman filter, and improve the performance of detecting faults through the residual chi-squared test [19]. However, the heavy calculation burden of artificial intelligence on the navigation computers limits its application in integrated navigation systems.
The variance shift outlier model (VSOM) combined with an extended Kalman filter (EKF) in this paper is utilized for the detection of GNSS measurement faults and the estimation of the variance shift for each measurement, which down-weights the measurement if required. The innovation of the EKF was fitted to the VSOM used to detect and identify the faulty measurements. The VSOM can be considered as the middle ground between FDI/FDE/FDR and robust estimation [20]. Unlike the methods of FDI/FDE/FDR [11,12], the size of the variance shift in the measurement determines the partial exclusion or inclusion of the measurement in the estimation in place of the complete deletion. Compared with robust filters [16,17], the identified faulty measurements are down-weighted, but not of all the measurements are weighted. The likelihood ratio (LR) test and score test statistics were used for the detection and identification of the faulty measurements, and the parametric bootstrap procedure was implemented for the significance assessment of both the LR and score testing and for multiple testing of the statistics.
The remainder of the paper is organized as follows. In Section 2, the mathematical models of the tightly coupled, INS/GNSS-integrated navigation system and variance shift outlier model are given. The novel method based on innovation and the VSOM are proposed in Section 3. Field test results, including the static test results and dynamic test results, are shown in Section 4. The summary for this paper is presented in Section 5.

Mathematical Model of a Tightly Coupled, INS/GNSS-Integrated Navigation System
A tightly coupled, INS/GNSS-integrated navigation system with a closed-loop errorstate, extended Kalman filter is described in this paper. The tightly coupled approach processes the raw pseudo-range GNSS measurement in a single filter, and the estimated position errors are utilized to correct the INS's navigation solution [21]. The architecture is shown in Figure 1. The linear dynamic equation for the error state in a time-varying linear system can be given as follows: The state vector in the filter is given by where , δv and δr are the error states of the attitude, velocity and position of the INS, respectively; δω b ib and δ f b are the gyroscope and accelerometer bias vectors; δC u and δĊ u denote, respectively, the GNSS receiver clock bias and the clock drift vector.
The linear observation equation is given by Assuming the discrete-time process, the system error dynamic equation and the observation equation can be rewritten as follows: where Φ k|k−1 , H k and z k are the transition matrix, the design matrix and the observation vector of the filter, respectively; ω k is the process noise vector, which is zero-mean Gaussian white noise with a process covariance matrix Q k ; and η k is measurement noise with covariance matrix R k , which is also zero-mean Gaussian white noise. The extended Kalman filter(EKF) algorithm consists of state updating and measurement updating. The prediction of statex k and its covariance matrixP k can be obtained by a state update as follows: When the observations are available,x k is the estimation of the state vector andP k is its covariance matrix. The measurement update in the filter is given bŷ where K k is the Kalman gain matrix, r k is the innovation vector and I denotes the identity matrix. The Kalman gain matrix is: The innovation vector is When a closed loop EKF is deployed, the value of the estimated state vector feeds back to the system, and the predicted state vectorx k becomes zero [12,21]. The innovation r k derives from the difference between the corrected pseudo-range vector ρ GNSS and the predicted pseudo-range vector derived from the solution of INS ρ I NS .
Under the fault-free condition, the innovation vector r k is Gaussian white noise with a zero-mean, and its covariance matrix P rk is given as: The horizontal alert limit (HAL) and the horizontal protection level (HPL) need to be compared with each other. When some faults in the measurements are not detected and identified, HAL and HPL could give the system protection. HAL is the maximal tolerable value of horizontal position error. When it is exceeded, an alert will be raised. The value of HAL is specified as 40 m in this paper. HPL is an upper bound. If the horizontal position error exceeds HPL, it shall be detected with a 99.9% probability. HPL needs to be computed in real time to check the integrity available and analyze the position-domain to determine whether the estimation could be used for the solution. Two computational passes are required for HPL in this tightly coupled, INS/GNSS-integrated navigation system, the first for checking integrity and the second for checking the availability of the final navigation solutions. A detailed calculation process for HPL is given in Section 3.4.

Variance Shift Outlier Model
The variance shift outlier model defined in this paper is the same as that in [20,22,23]. The linear model is as follows: for x is a q × 1 parameter vector and η is a p × 1 random error vector that is assumed to obey a Gaussian distribution with a zero mean and the variance σ 2 I. The residual maximum likelihood (REML) ignoring constants LR takes the form wherex = (H T H) −1 H T z is the best linear unbiased estimate (BLUE) of x under the linear model (14). The REML estimateσ 2 = (z − Hx) T (z − Hx)/(p − q) of σ 2 is unbiased. Suppose the ith measurement has an inflated error variance. i is the number less than or equal to the dimensions of z, which denotes the position of the fault in the measurements. The measurement has error variance σ 2 ( i + 1); i ≥ 0, and i denotes the inflated factor in the variance of the ith measurement. A variance shift model for the ith measurement takes the form where d i is the ith p × 1 unit vector with 1 in the ith element and 0 elsewhere; δ i is a random coefficient with zero-mean and its variance is i σ 2 where i ≥ 0. W i = i d i d T i + I is a diagonal matrix with i + 1 corresponding to z i and 1 elsewhere. I is a p × p identity matrix. The variance matrix of z the data under a VSOM model (16) is The variance corresponding to the ith element of z inflates i + 1 more than the variances of the fault-free elements.
The REML log-likelihood function (REML LLF) for the ith VSOM (16) can be expressed as follows: where . The REML LLF for the ith VSOM can be given in another way: where . The REML estimates of variance and inflated factor are [24] whereσ 2 0 =η Tη /(p − q) and t 2 i = s 2 i /σ 2 0 are the error variance estimate and the Studentized residual with degrees of freedom p − q for the ith measurement under the null model (14), respectively. Since t 2 i /(p − q) is distributed as beta(1/2, (p − q)/2), t 2 i should less than p − q. In the case of t 2 i = p − q, the estimate ofˆ i would approach infinity [25].

Methodology
A linear model that follows is used here.
where z k is the innovation of the Kalman filter, which is the same as the one defined in Equation (12). x f rame is the vector of position errors and receiver clock biases, whose dimensions are equal to three plus the number of observed constellations. G k is the geometry matrix defined in [21,26].

Establishment of Test Statistics
Normalize the measurement equation as follows (22) [21,27]: where: P −1/2 rk is inverse of the upper triangular matrix by the Cholesky decomposition of P rk . Then Equation (23) has the same form as Equation (14). According to Equations (14)- (21), some parameters under fault-free conditions can be rewritten as where N Sat is the number of observed navigation satellites which is equal to the dimensions of z k ; N Sta = N Const + 3 is the dimensions of unknown parameter x f rame ; and N Const is the number of observed constellations.η i contains the inflated error variance of the ith satellite measurement; see Appendix A. To confirm whether there is a variance shift parameter i larger than zero in the ith measurement, testing of the hypotheses is established as follows: The alternative model is fitted under the VSOM of the ith measurement.

Likelihood Ratio Test Statistics
To test the hypothesis, the likelihood ratio test (LRT) statistic is established by [20,23] whereψ 0 = (0,σ 2 0 ) T is a variance estimate under the null hypothesis andψ i = (ˆ i ,σ 2 0 ) is one under the VSOM alternative hypothesis of the ith measurement. RL 0 (ψ 0 ; z nk ) is the REML LLF for z nk atψ 0 , and RL i (ψ i ; z nk ) for z nk atψ i . The subscript i denotes the ith VSOM.
Under the null hypothesis, the REML LLF could be given by and under the ith VSOM alternative hypothesis, the REML log-likelihood function is given by Hence, LRT is a monotonically increasing function of t 2 i .

Score Test Statistics
There is an asymptotical equivalence between the LRT and the score test under the null hypothesis. However, there is a computational advantage to the score test statistic. The score test statistic does not need to fit the model specified under the alternative hypothesis, and it is only required to estimate under the null model [20,22]. A score test statistic calculation only requires the score vector and the information matrix under the null hypothesis. For the variance parameter i , the score function by differentiating Equation (19) with respect to i is given by Under the null hypothesis, i = 0 and σ 2 =σ 2 , the score function can be expressed as follows.
The negative expected value of the second derivatives of (19) is specified as the information matrix for i under the ith VSOM, and the information matrix can be partitioned as follows: Using Equations (23)-(26), the expected information matrix under the null hypothesis is given by For testing H 0 : i = 0 against H A : i > 0, the score test statistic takes the form The score test statistic can be rewritten as

Significance and Multiple Testing
There are two situations in the practical scenario [23]: a measurement is suspicious and identified as a possible fault before analysis; the measurements are screened without prior information about the faulty measurements. In the second situation, ideally, all potential anomalous measurements are identified and included into the modeling process.
To determine how many and which measurements are faulty, the LR and score test statistics for each measurement are calculated. Then, testing more than one hypothesis simultaneously becomes the primary problem. For the appropriate sampling distribution of the LR and score test statistics, a parametric bootstrap method is implemented for significance assessments of the tests and multiple testing of the statistics. The parametric bootstrap procedure is can be divided into five steps: S1: Estimate (23) under the null hypothesis, obtaining parameter estimatesx andσ 2 .

S2: Generate new measurements
where η * nk is simulated as N(0,σ 2 I n ). S3: Fit the null hypothesis (23) and obtain bootstrap LR test statistic LRT * i ; i = 1, . . . , N Sat and score test statistics S * i ; i = 1, . . . , N Sat . Obtain the order statistics from each set. S4: Step 2 and step 3 are required to be repeated B times, for B is reasonably large. An empirical distribution(ED) of size B for each order statistic is generated. S5: Calculate the 100(1 − α)th percentile for each order statistic for the required α, where α is the significance level. When the LR test statistics exceed their respective thresholds (the score test is the asymptotic equivalent to LRT), these measurements are deemed to be all faulty and should fit the modified model, including a variance shift. The measurements that are deemed to be faulty will be down-weighted in the estimation procedure.
The flow chart of detecting faulty measurements based on the VSOM approach is shown in Figure 2, where D is a N Sat × N F matrix that composes N F vectors with 1 in the ith position and 0 elsewhere, δ is a N F × 1 vector with random elements and N F is the number of faulty measurements identified.

Down-Weighting
After obtaining the inflated factor of the variance shift in the measurement according to Equation (21), the down-weighting process is implemented as follows: where the new covariance matrix P n rk of the innovation is obtained by multiplying the inflated factor by P rk , and P rk is defined in (13). P 1/2 rk is the same as when defined in (24). W i = I + D T ΠD; Π is a diagonal matrix with i in the elements corresponding to the faults, and 0 elsewhere. Then, substitute P n rk into the Kalman filter in place of P rk .

Horizontal Protection Level Computation
The computation of HPL is described as follows [28,29]. HPL1 was set to 5.33σ H , where σ H was obtained from the portions of covariance matrix that corresponded to the state error vector elements about horizontal position errors, and 5.33 was derived from the missed detection rate P MD of 10 −3 /h. The calculation of HPL2 is similar to that of traditional receiver autonomous integrity monitoring (RAIM) HPL [30]. HPL2 obtains the protection level by projecting the test statistic to the position domain. The ratio is referred to as SLOPE. The calculation of HSLOPE is given as follows: where K n 7i and K n 8i are the gain coefficients of biases in the measurement derived from gain matrix K k corresponding to the north and east directions, respectively. S ii is the ith element of diagonal matrix (I − H k K k ) T (I − H k K k ).
HPL2 can be calculated by where P bias is the square root of the non-centrality parameter corresponding to a missed detection rate equal to 10 −3 /h. HPL can be given by

Field Test
In this section, the proposed method is proven to be feasible by the a static test and a vehicular dynamic test of the tightly coupled, INS/GNSS-integrated navigation system. The fiber-optic-gyroscope-based SINS in both two experiments was independently developed by Harbin Engineering University. The sampling rate for the IMU was set at 100 Hz. The parameters of the INS are shown in Table 1. The GNSS receiver used in both two experiments was Unicore UB370. The rising edge of 1PPS was used as the trigger signal to reset the counter in the navigation computer. The counter give the raw date from IMU time information in the interval between two trigger signals. The sampling frequency of IMU was high enough, and the offset of the local system clock was not greater than 0.01 s. The output rate of receiver was set at 1 Hz. The number of bootstrap samples B was set at 1000. The significance level α was set at 0.01 [31].

Static Test
The INS was placed on a marble isolation table, and the antenna of the receiver was placed outdoors. The GNSS receiver tracked dual-constellation navigation signals: GPS L1: 1575.42 MHz; BDS B1: 1562.098 MHz. The raw navigation data from the IMU and GNSS receiver were collected for 21,600 s. The arrangement is shown in Figure 3.  Table 2.  Table 2 shows that the ith measurement is just a potential fault when t 2 i > 1. The faults were confirmed after either an LR test or a score test. Only the eighth measurement was the fault, and the other measurements were fault-free.
The comparison of the position solutions with and without fault detection and identification based on VSOM is shown in Figure 4. The timespan of these data was six hours. There were at least four errors in the range of the 500th second to the 2000th second, obviously before fault detection and identification, and the errors caused an abrupt jump over 5 m. After the detection and identification process, the errors disappeared. The RMSE of the position before fault detection and identification was 5.8459 m, and that after fault detection and identification was 4.4656 m. Considering the position errors before fault detection and identification, this indicates that the proposed method based on the VSOM can detect and identify the faults correctly, and the faulty measurements in the estimation ensure the precision and reliability. The proposed method exhibited superior performance for the detection and identification of faults during static testing.

Dynamic Test
A vehicular dynamic test was performed in an urban environment. The test equipment is shown in Figure 5. The IMU and the receiver antenna were fixed on the vehicle. The GNSS receiver tracked dual-constellation navigation signals: GPS L1: 1575.42 MHz; BDS B1: 1562.098 MHz. Both GPS and BDS were employed in this experiment, and the number of observed satellites was the sum of that of two systems. Figure 6 shows the positional dilution of precision (PDOP) and the number of satellites visible during the test. The average of the PDOP was 2.02, and the average number of observed satellites was 10.67.  Figure 7 shows an overview of vehicle trajectories, and shows the difference between the trajectories before and after fault detection and identification based on the VSOM. When one or more faults occurred, the trajectory became deformed. As shown in the figure below, a fault occurred at the 350th second, and the fault caused an abrupt disturbance. The proposed method was implemented, and the trajectory was not affected by the fault. Obviously, after down-weighting the fault in the estimation process, the trajectory maintained stability. To further reveal the details, the distance from the initial position to the current position over time is shown in Figure 8. As shown in the figure, when the faults occurred, the east/north/upward directions all experienced some disturbance. Before fault detection and identification, the position moved 10 m to the east and 5 m to the south abruptly. A jump in the upward direction of 20 m occurred. After the fault detection and identification, the fault was down-weighted in the estimation process, and the solution was not disturbed. The transitions from the previous position to the fault position, and the fault position to the next position were seamless. This result indicates that the fault detection and identification based on the VSOM could eliminate the vibrations in the estimation caused by the faults. When the weight factors were used for down-weighting approaches to infinity, the corresponding measurements in the estimation were nearly deleted.

Conclusions
A fault detection and identification method for tightly coupled, INS/GNSS-integrated navigation systems is proposed in this paper. The method emphasizes the faults in GNSS measurements. This method is beneficial to ensuring the precision and reliability of a tightly coupled, integrated navigation system. The VSOM employed in this paper can be viewed as a trade-off between a faulty detection and exclusion algorithm and robust estimation. The size of the variance shift in the measurement is the key to partial exclusion or inclusion of the measurements in the estimation process. Either the LR statistics or score test statistics are presented as objective measures for determining whether the measurements are faulty. Parametric bootstrapping was employed in this paper to assess the significance and the handling of multiple testing. The performance of the proposed method was verified by a static field test and a dynamic field test. As shown in the results, the faults occurring during the navigation process can be detected and identified accurately. After down-weighting the fault with proper weight factors in the estimation procedure, the performance and precision can be ensured. However, the fault detection and identification based on the variance shift outlier model still need to be improved. More applicable scene tests and more flexible test criteria will be explored to perfect the method in the future.