Analysis of Single and Double Faults Direction and Magnitude in Measurement and State Models of Tight GPS/INS System

Analysis of Single and Double Faults Direction and Magnitude in Measurement and State Models of Tight GPS/INS System

DOI: https://doi.org/10.52939/ijg.v19i7.2745

Integration of Global Navigation Satellite Systems (GNSS) with Inertial Navigation Systems (INS) is widely utilized as the best navigation solution in various applications such as military and civil aircraft, aerial photogrammetry, and Mobile Mapping Systems (MMS) [1]. Such integration is necessary to overcome the shortcomings of a standalone GNSS or INS System because both systems have complementary error characteristics. Therefore, a GNSS/INS integrated navigation system offers improvements in the availability and reliability of the full navigation solutions [2]. As GNSS/INS integration is usually performed using an Extended Kalman Filter (EKF), the outcomes of EKF may include large errors due to different faults in either GNSS or IMU system. Due to different sources, these faults have different characterizations in terms of magnitude, direction, and occurrence (prompt or slow growing). Over the last two decades, several Fault Detection and Exclusion/Identifications (DE/FDI) schemes have been developed to localize and exclude faults either in measurement or system state models of the GNSS/INS system. These schemes evolved from hypotheses testing for deciding the accepted region where a fault does not influence the navigation solution. Otherwise, the rejected region is decided if the size of the fault exceeds the predefined critical value. Receiver Autonomous Integrity Monitoring method RAIM which is used to check the reliability and integrity of the GNSS/INS system [3], Autonomous Integrity Monitoring by Extrapolation (AIME) [4], Innovation-based, and residuals-Based methods [5], Solution Separation (SS) method [6] and [7], Quality Control (QC) [8], Extended RAIM (ERAIM) [9] are examples of these schemes.

Faults in measurement and system state models of the GNSS/INS system have distinctive characteristics in terms of size, direction, and sudden or slowly growing faults. The size of the faults can be represented as a scaler and obtained from fault test statistics whereas the direction is an angle between two fault vectors [10] and [11]. Additionally, prompt faults may occur in IMU dynamic model due to sudden changes in vehicle trajectory while ramp faults are tiny fault that cumulates in GNSS pseudorange observations or system states with time propagation. Therefore, investigating the behavior of the faults is necessary to detect, identify, and remove them using a suitable algorithm that guarantees continuity and availability throughout the navigation time.

Intensive research works e.g., [2] [9] [12] [13] and [14] have already been performed to detect, identify, and exclude faults in integrated GNSS/INS systems. To do that, the detection, identification, and adaptation method (DIA) based on a single epoch was used. In this method, a global test is initially used to detect the faulty epoch and if an abnormality is detected in a certain epoch, the local test is conducted to identify the faulty measurement/measurements in that epoch. Then, the faulty measurements are removed in the last step (adaptation step). Another research group e.g., [15] [16] [17] [18] and [19] assumes that any fault that occurs in the past continuously influences the current-time innovations and state estimation of the extended Kalman filter (EKF). An example of these methods is the autonomous integrity monitoring by extrapolation (AIME) method [15]. The AIME methods are based on chi-square detection in which the measurements used are not limited to a single epoch.

Identifying the actual faulty measurements and excluding them is a crucial step because wrong identification may cause discontinuity/unavailability of the navigation solution. On some occasions, wrong identification occurs due to the high correlation between the actual faulty measurement and other measurements. Thus, the correlation between fault test statistics is deemed as a measure of separability as such a high correlation means a low probability of fault separability and vice versa. Correlation between fault test statistics when the test is executed to identify a single fault at a time has been investigated [9] [10] and [11]. For multiple simultaneous faults, the multi-dimensional correlation was pioneered by Fortner [10] and then developed by Li [20]. Wang et al., [21] and Almagbile [11] investigated multi-dimensional correlation for multiple faults in GNSS/INS system. However, more detailed analysis is required to understand the probability separability when multiple faults occur in GNSS/INS measurement and/or dynamic models. In this research, the behavior of multiple faults in pseudorange measurement and predicted states of tight GPS/INS system is investigated. More specifically, this research aims to (i) analyze the magnitude and direction of multiple faults in both measurement and dynamic models, and (ii) develop correlation coefficients to analyze the cases of faults separability. The rest of this paper is organized as follows. Section 2 describes the snapshot approach for fault detection and identification in GPS/INS system. This includes hypotheses testing, Identifiability, and separability with multiple directions and magnitude. Section 3 provides developed correlation coefficients between fault test statistics. Simulation, results, and analysis are presented in Sections 4 and 5 respectively. Section 6 includes the concluding remarks.

2. Snapshot Approach for Fault Detection in GPS/INS system

For tight GPS/INS integration, an extended Kalman filter (EKF) is normally used. After the linearization of EKF, the dynamic and measurement model of the linear discrete-time system can be described [2] and [22]:

Equation 1

Equation 2

Where:

x_k is state vectors at time k

x_k-1 ∈ ℜ ^n×1 is state vectors at time k-1

Φ_k-1 ∈ ℜ ^n×n is transition matrix

z_k∈ ℜ ^m×1 observation vector at time k

Z_k ∈ ℜ ^m×n is observation matrix

w_k-1 ∈ ℜ ^n×1 is uncorrelated white process

v_k ∈ ℜ ^m×1 is measurement noise with zero mean with covariance matrices Q_k and R_k

To utilize the snapshot approach for quality control in an integrated GPS/INS system, the Kalman filter is converted into least squares. The least squares navigation solution consists of the measurements and the predicted states of the Kalman filter and can be written as [23] [24] and [25]:

Equation 3

Equation 4

The observation vector l_k is formed by combining the observation Z_k and predicted states X̄_k vectors. The design matrix A_k is formed by the measurement matrix H_k and identity matrix I_(n×n), v_k is the residuals vector contains the residuals vectors of measurements v_zk and predicted states v_yk . The variance-covariance matrix D_k of the measurement vector l_k can be written as [21]:

Equation 5

Where D_k is formed from the measurement noise covariance matrix R_k and the predicted states covariance matrix P̄_k of Kalman filtering, P_k^-1 is the weight matrix. The optimal estimate of the state parameters x̂_k and error covariance matrix Q_x̂k can be written as [14]:

Equation 6

Equation 7

The residuals vector v_k and cofactor Q_vk can be written as:

Equation 8

Equation 9

The least squares approach can be then used for fault detection and identification of the integrated GPS/INS system. For simplification, the symbol k which represents the epoch number is omitted from the subsequent equations.

3. Detection, Identification and Separability of Multiple Faults

3.1 Fault Detection

In this step, potential faults are detected using the global model test to check whether the system includes faults or not. To do that, hypothesis testing is carried out by testing a null hypothesis against an alternative hypothesis as follows [12] and [26]:

Equation 10

Where:

H₀and H₁ are the null and alternative hypotheses respectively, and are the posteriori and priori variance respectively. nand u are the known and unknown parameters where the difference between them is known as the degree of freedom (f).

The global model test is then can be computed as:

Equation 11

Where χ²is the Chi-Square distribution, α(gf) is the critical value of the global model test based on a significant level and degree of freedom 1-α.

3.2 Fault Identification

In the case that the global model test detects an abnormality in a certain epoch, a local test is employed to identify the measurements that hold faults. Once again, the null and alternative hypotheses are tested [10]:

Equation 12

Equation 13

Where ∇_i is the estimated size of an outlier in the i ^th observation, h_i=[0, …0, 1, 0] ^Tis a fault vector in which the i ^th element is equal to one and all other elements are equal to zero. Assuming that there is an outlier in the i ^thobservation, the alternative hypothesis is accepted and the residual vector can be written as:

Equation 14

To identify a fault in the i^th observation, a test statistic is implemented based on a data snooping approach which means testing each observation individually to identify the observation that holds a fault. This test is written as [27] and [29]:

Equation 15

Where w_i is the test statistics, q∇_i=(h^T_iPQ_vPh₁)^-1 . In the case that the null hypothesis is accepted, has a standard normal distribution. However, if the alternative hypothesis is accepted, this means that there is a fault in i ^th observation and is thus has a non-centrality as follows:

Equation 16

For one tail test, the critical value can be determined from the significant level 1-α as:

Equation 17

With a given power of the test (1-β) and false alarm (1-α), the minimal detection bias (MDB)∇δ_i can be computed as [13] and [26]:

Equation 18

Where δ₀ is a non-centrality parameter

For multiple faults, the test in Equation (15) is generalized [11] [12] and [26] as:

Equation 19

Where θ is the number of faults H_i is the fault matrix with rank θ, containing zeros with a one in each column corresponding to the number of faults in the observation vector, 1-α_Tiθ is the level of significance for multiple faults.If on some occasions the alternative hypothesis is accepted, T^θ_i become non-central normal distribution and then the non-centrality parameter λ can be written as [10] [28] and [29]:

Equation 20

Where ∇δ is the true fault vector. The minimal detection bias (MDB) in this case can be computed by determining the multiple correlation coefficients. The multiple correlation coefficients are equivalent to the correlation coefficients between two single-fault test statistics. The multiple correlation coefficients ρ^θ_ij and MDB ∇₀δ^θ_i can then be written respectively as [11] and [26]:

Equation 21

Equation 22

with bounds of ρ^θ_ij are: 0 ≤ ρ^θ_ij ≤ 1

3.3 Faults Separability

When one rejects the null hypothesis in favor of the alternative hypothesis, this means that the test statistic is larger than the critical value. This is normal when single or multiple faults occur in the observation. However, a wrong decision may happen when one rejects a proper alternative hypothesis and accepts the other and this lead to type III error [10]. The reason behind accepting the wrong alternative hypothesis and rejecting the true one is due to the correlation between fault tests. This correlation can be considered as a measure of separability as when the correlation is high, the probability of correctly identifying the true fault position is low, and vice versa. The correlation ρ_ij between w_i and w_jcan be written as [10] and [30]:

Equation 23

Likewise, the correlation between T_i and T_j can be generalized from the nonsquare matrix [10] and [21]:

Equation 24

From Equatiom (24), the correlation coefficients between the two tests can be expanded as:

Equation 25

Where:

Equation 26

Equation 27

Equation 28

Equation 29

When double outliers test is conducted, the matrix M₁₂ is a quad matrix that includes multi-dimensional correlation between test T_i and T_j. Now the correlation is equal to the maximum eigenvalue λ ofM₁₂as [10]:

Equation 30

The correlation presented in Equation (26) can be rewritten with the geometric form [10] and [11]:

Equation 31

Where ε is the angle between two fault vectors, and it can be calculated as:

Equation 32

The angle ε behaves oppositely to the ρ²₁₂ as such a high correlation between tests reflects a small angle and vice versa.

4. Simulation and Experimental Scenarios

4.1 Simulation of Tight GPS/INS System

In this research, tightly coupled GPS/INS integration was simulated using GPSoft ® Toolbox software. The simulation started by creating a reference trajectory with zero (0,0,0) initial values of position, velocity, and attitude. The ground truth of INS data was created without errors. Errors in gyros, accelerometer, initial alignment, and velocity were then added to the INS data (Table 1). GPS pseudo-range observations were derived from five visible satellites during the navigation time (649 seconds) (Figure 1).

Table 1: Simulated INS error sources

Figure 1: Skyplot of visible satellites

The GPS pseudo-range and INS data with and without errors were integrated using an extended Kalman filter. The dynamic model of the Kalman filter consists of 17 states:

Equation 33

Where:

r_E, r_N and r_U are east, north, and up position, respectively

v_E, v_N and v_U are velocity in east north and up, respectively

φ_E, φ_N and φ_U are altitude in east, north, and up, respectively

ε_x, ε_yand ε_z are gyros errors X, Y, and Z, respectively

∇_x, ∇_y and ∇_z are the accelerometer errors X, Y, and Z, respectively

X̂_b is receiver clock bias

X̂_f is receiver clock drift

The INS output rate was 200 Hz and it was updated every 1 Hz with pseudo-range observations. To show the performance of the integrated system, position errors in East, North, and Up components were computed from the difference between the GPS/INS integration and the ground truth data. The INS-unaided, and INS/GPS relative to the ground truth trajectory are shown in Figure 2. The magnifier window illustrates the difference between the INS-unaided, GPS/INS, and ground truth vehicle trajectories in detailed scale.

4.2 Faults Simulation and Experimental Scenarios

After performing the GPS/INS integration, faults with different magnitudes were injected in pseudorange observations as well as in XYZ position predicted states. These simulated faults are presented in Table 2. The experiments were conducted to illustrate the magnitude and direction of faults using these scenarios:

· A fault of 15m was injected in a single pseudornage observation and X position state. Likewise, a fault of -15m was injected in a single pseudornage observation and X position state.

· Two faults of 15m were injected in two pseudornage observations and X and Y position states. Then two faults of -15m were injected in the same pseudornage observations and position states.

5. Results and Discussion

5.1 Single Fault Test in The Measurement Model

Single fault test along with single correlation coefficients are depicted in Figure 3(a) and Figure 3(b). From Figure 3(a), one can notice that when a fault of 15m was injected in pseudorange 1 the fault test values identified this measurement as faulty measurement. Also, the influence of the outlier on other measurements can be noticed particularly in measurement 3. Unexpectedly, the fault test value of measurement 3 when a fault occurs in measurement 1 becomes less than those without an outlier. This explains that some measurements do not necessarily follow the same fault direction. It can also be noticed that the influence of the fault on measurements 2, 4, and 5 is tiny. The correlation between the fault test of measurement 1 and all other measurements is attached to Figure 3(a).

Figure 2: 2D and 3D vehicle trajectories include INS-unaided, GPS/INS relative to the ground truth trajectory

Table 2: Simulated faults scenarios in pseudorange and predicted states

Figure 3: Single fault of +15m and correlation analysis for (a) pseudorange 1 and (b) pseudorange 2

Figure 4: Single fault of -15m and correlation analysis for (a) pseudorange 1 and (b) pseudorange 2

As can be seen, all correlation values are less than 0.2 and this means that the influence of the faulty measurement on other measurements is tiny. In Figure 3(b), a fault was injected in measurement 2 and it is clear that the fault test value of this measurement is the highest among others but it is still less than the critical value (10.82) [13]. This means that the value of the fault test varies in accordance with the precision of the measurements.

Figure 4(a) shows a single fault test when a fault of -15m is injected in measurement 1. As can be seen, measurement 1 is identified as a faulty measurement but its value is less than those in the previous fault case (+15m). It can also be noticed that all measurements are influenced by the fault and almost all of them have a fault value less than those without fault. In Figure 4(b), when a fault of -15m was injected in measurement 2, the situation becomes different from those in Figure 3(b) because the test value exceeds the critical value.

5.2 Single Fault Test in the Dynamic Model

In this scenario, a fault of 15m and -15m is simulated in the X position state and Y position state as can be seen in Figures 5(a), (b), (c) and (d) respectively. When the magnitude of the fault is 15m, the test can identify that the X position (Figure 5(a)) and Y position (Figure 5(c)) are faulty states. In Figure 5(b), a fault of -15m occurs in the X position and the test reflects a large value. Additionally, both the X position and Y position influence each other more than the Z position state. This is due to the higher correlation coefficients between the tests for the X and Y states compared to that between the X and Z states and also between the Y and Z states (Figure 6(a), (b) and (c)).

5.3. Double Fault Test in the Measurement Model

In this scenario, a fault test is implemented under the presence of two faults simultaneously. Figure 7 illustrates the situation of double faults in the combined pseudorange (1 and 2) and compared it when these measurements do not have any faults. As can be seen, the fault test value of the measurement pair (1,2) before and after 15m faults was injected in both measurements is approximately 1.5 and 34 respectively. This confirms that the fault test was successful in identifying the faulty pair. Additionally, it can be seen that the combinations (1,3), (1,4), and (1,5) exceed the critical value (13.82) too [13]. This was because those combinations include one faulty measurement, and thus, the faulty measurement adversely affects all the measurements that are combined with it. Nevertheless, despite that the combinations (2,3), (2,4), and (2,5) include a faulty measurement, they pass the test because their values are less than the critical value.

Figure 5: Single fault of 15m and -15m in X, Y, and Z position states respectively

Figure 6: Single correlation analysis for (a) X position, (b) Y position and (c) Z position states

Figure 7: Double faults of +15m (a) and -15m (b) in combined pseudorange 1 and pseudorange 2

Figure 8: Double fault of 15m (a) and -15m (b) in combined X and Y position states

This situation is found earlier in this research when a single fault does not identify measurement 2 as a faulty measurement even when a fault of 15m was injected in this measurement. The test for the pair (3,4) becomes less than that without faults in the pair (1,2). When a fault of -15m was simulated in the pair (1, 2), the fault test identifies those pairs as faulty pairs but with higher values than that of the +15m fault case. Notably, the values of all the measurements which are combined with measurement 1 become less than in the previous case. Furthermore, it can be seen that the fault test for the pairs (2, 3), (3, 4), (3, 5), and (4,5) reveals the opposite situation to the +15m case.

5.4. Double Fault Test in The Dynamic Model

Similar to the previous scenario, two cases of faults were injected in X and Y position states. From Figure 8, it can be seen that all combined states (X, Y), (X, Z), and (Y, Z) exceed the critical value. Additionally, in the case of a -15m fault, all the combinations values are higher than that of a +15m fault case.

5.5. Fault Direction and Multi-Dimensional Correlation

5.5.1 Measurement model

For a better understanding of the influence of the measurement that holds fault on other pure measurements, fault test vectors before and after injecting a fault of +15m in measurement 1 along with the correlation coefficients between fault tests for each measurement are shown in Figure 9. It is clear that when a fault occurs in measurement 1, the direction of the fault vector becomes different from that without fault. The faulty measurement pulls each measurement in a different direction and this depends on the correlation between the faulty measurement and other measurements. For instance, while measurement 2 is pulled in the same direction as the faulty measurement, measurement 3 is pushed in the opposite direction. The direction which is represented by an angle between the x-axis and the correlation coefficient vector is shown in Figure 9(c). The higher the correlation between fault tests is, the smaller the angle becomes and vice versa.

Figure 9: (a) Fault direction without fault, (b) with a single fault of +15m in SV1 and (c) the corresponding angle of the correlation vector between fault test

Thus, the maximum correlation between the fault tests for the measurements 2 and 4 reflects the smallest angle. The correlation vector between measurements 1 and 2 as well as measurements 1 and 4 did not appear because they are hidden behind the correlation vector between measurements 2 and 3. Figure 10 shows the situation when two faults of +15m occur in measurements 1 and 2. From the Figure it can be seen that the combination (1, 2) pulled all the measurements that have common observations with either 1 or 2. It is clear that the combinations (1, 3), (1, 4), (2, 3), and (2, 4) have higher values compared with that of the (3, 4) combination.

To clarify the relationships between fault tests for the measurement combinations, the correlation coefficient, and the corresponding angle under three fault cases are depicted in Table 3. The correlation coefficient includes a quad matrix (M) (refer to Section 3.3) which represents the interaction between two parts:

The general relationship between the combinations indicates the magnitude and direction of any combination related to the faulty ones.

The specific relationship between the combination elements itself and also between the elements of different combinations.

Three correlation cases can be seen as follows:

Case 1: The correlation between the combination (1, 2) and itself is always equal to 1. However, when this combination was divided into four relationships as the relation between (1 and 1), (2 and 2), (1 and 2), and (2 and 1), their correlation values are 1, 1, 0.102, and 0.102 respectively. Thus, the corresponding angles are equal to 0°, 0°, 90°, and 84.14°. Notably, the angle is 90° and 0° when the correlation is 0 and 1 respectively.

Case 2: The two fault vectors include one common measurement. As the combination (1, 2) and (2, 3), the measurement 2 is common between those combinations and thus the correlation is 1 with a corresponding angle equal to 90°. It can be noticed that the correlation between (1 and 2), (1 and 3), (2 and 2), and (2 and 3) is 0.102, 0.063, 1, and 0.065 respectively. The corresponding angles are 84.15°, 86.42°, 0°, and 0, 86° respectively.

Figure 10: (a) Fault direction without fault and (b) double faults of +15m in SV1 and SV2

Table 3: Correlation and the corresponding angle for individual and combined measurements

Case 3: when no common measurements occur between two combinations. For example, the combination (1,2) and (3,4) have no common elements and thus the correlation and the corresponding angle are 0.126 and 82.74° respectively. Tiny correlation values (ranges between 0.06 and 0.09) can be seen between (1 and 3), (1 and 4), (2 and 3), and (3 and 4) with high corresponding angles range between 84.51° and 86.41°. The direction of the fault vectors is visualized in Figure 12.

5.5.2 Dynamic model

Figure 11 shows fault test vectors before and after injecting a fault of +15m in the X position state along with the corresponding angles of the correlation coefficients between fault tests. When a fault occurs in the X position state, the direction of the fault vector for all states (X, Y, Z position, and X velocity states) behave differently from that without fault. What is more, the correlation plays a crucial role in determining the magnitude and direction of the fault vectors. The faulty state (X position) pulled the Y position state to its direction because of the relatively high correlation between those states. On the contrary, the Z position and X velocity states have a tiny correlation with the faulty state and thus their directions and magnitude are different from that of the Y position state. The maximum correlation coefficient is equal to one as shown in the Figure 11.

When two faults of +15m occur in the X and Y position states as shown in Figure 12 the faulty states pulled the X and Z position state and the X position and the X velocity states to their direction and magnitude. Note that the fault vector of the combination (X, Z) position states is hidden behind the vector of (X position, X velocity) states because both combinations have very close fault magnitude.

Table 4 illustrates the nested correlation coefficient and the corresponding angle under three fault cases. As shown previously, it can be seen that the correlation between the states (X, Y) and itself is always 1 with an angle equal to 90°. However, the individual correlation between X and Y is around 0.36 with an angle of 68.74°. From Table 4 and Figure 12, the correlation of the combination (X, Y) and (X, Z) is equal to 1 (maximum correlation) where the correlation between X and Y is 0.36, X and Z (0.11), Y and Y (1), Y and Z (0.02) with corresponding angle 68.74°, 83.52°, 0°, and 88.77° respectively.

Figure 11: (a) Fault direction without fault, (b) with a single fault of +15m the X position state and (c) the corresponding angle of the correlation vector between fault test

Figure 12: (a) Fault direction without fault and (b) double faults of +15m in in X and Y position states

Table 4: Correlation and the corresponding angle for single and combined states