Abstract

The problem of fault detection and isolation (FDI) on inertial measurement units (IMUs) has received great attention in the last years, mainly with growing use of IMU strapdown platforms using fiber optic gyros (FOG) or micro electro mechanical systems (MEMSs). A way to solve this problem makes use of sensor redundancy and parity vector (PV) analysis. However, the actual sensor outputs can include some anomalies, as impulsive noise which can be associated with the sensors itself or data acquisition process, committing the elementary threshold criteria as commonly used. Therefore, to overcome this problem, in this work, it is proposed an algorithm based on median filter (MF) for prefiltering and chi-square cumulative sum (-CUSUM) only for fault detection (FD) applied to an IMU composed by four FOGs.

1. Introduction

The design of a FDI algorithm for applications on IMUs can, in general, be divided in two types.The first one, named analytical redundancy, takes into account a mathematical model of the system in which the IMU is used. This method generates a residue vector as a result of the state observer [14] in order to indicate the IMU operational status. The second one, named sensor redundancy, is based on specific geometrical configurations using extra sensors and solving the FDI problem with the aid of the parity vector (PV) analysis [5, 6]. However, the actual sensor outputs can include some anomalies, as impulsive noise which can be associated with the sensors itself or with the data acquisition process [7], violating the elementary threshold criteria as commonly used. This work applies the second type of algorithm to the FD problem. The use of four sensors (minimal redundancy) constrains the algorithm to detect the fault not allowing the identification (isolation) of the faulty sensor [8]. To overcome the anomalies addressed before, in this work is proposed an algorithm based on MF that performs the prefiltering of FOG outputs and the use of the -CUSUM for FD problem.

2. Background

2.1. Geometry

The geometrical arrangement used in this work considers four gyros mounted on the faces of a tetrahedral structure (tetrad), and three accelerometers in a triad configuration internally fixed in the tetrahedral. The analysis performed here takes into account the gyros only, and the extension for accelerometers is straightforward. The tetradconfiguration and the reference frame are shown in Figure 1. The mathematical representation of the arrangement of the gyros is given in terms of direct cosine matrix (DCM), and is obtained from angular relationship between sensor axes and the analytical triorthogonal axes (analytical triad). Therefore, considering the schema shown in Figure 1, where , the DCM is given by The matrix relates the sensor measurements with the angular rate in the main axes in the form of The estimate of the angular rate components in the main axes can be obtained from (2.2) in the following manner: where is the vector of the gyro outputs, is the angular rate vector about main axes, is the angular rate estimate vector, and is the generalized inverse of .

Equation (2.4) provides the best state estimation in the least squares sense.

2.2. Parity Vector

The sensor equation considered in (2.2) can be rewritten by adding the faults, biases, and noises components as follows: where is the bias vector whose magnitude is constant, is the fault vector, and is a random term vector (Gaussian noise). Applying the singular value decomposition (SVD) on , it can be obtained the range and null spaces from this matrix [9]. In addition, also is computed the biases influence on the arrangement. Decomposing as follows: where , , and are matrices obtained from SVD of . The matrix is a diagonal matrix whose elements are the eigenvalues of . The superscript () indicates the transpose. Applying (2.7) into (2.5) and multiplying both sides by , it can be obtained the following relationship: Partitioning as follows: where and , and applying it into (2.8), the resulting equations are Equation (2.10) leads to least squares estimate of , what is equivalent to (2.3), and can be expressed by The meaning of (2.11) is that, if sensors faults and biases are zero, the resulting product of parity vector with sensor measurements is a white noise with zero mean. Otherwise, if the biases and/or faults values differ from zero, (2.11) is a “weighted’’ sensor errors summation. Then, (or ) is null space of and is the parity vector.

2.3. Median Filter

In the image processing field, it is very common to employ filters that preserve edges or abrupt transitions between distinct parts of the image or remove salt and pepper noise. These filters generally compare the pixel under observation with its neighbors at certain window and take a decision based on a statistical or threshold criteria. The simplest filter that meets these requirements is the median filter (MF). Being the MF not an optimal filter, it preserves discontinuities as jumps [10] and is also a robust estimator [11], besides being an easy and fast way to remove outliers from signals [7].

2.3.1. Recursive Median Filter

In the time sequence processing by MF, it is possible to process the samples in the filter window in two ways. (a) Process the samples in the window to obtain the output at discrete time (), shift the window to process the time () over the original data without replacement; this is the nonrecursive processing. (b) Process the samples in the window to obtain the output at discrete time (), replace the original sample by the filter output , shift the window and process the new samples with replacement; this is the recursive median filter. The nonrecursive MF is shown by (2.14) and the recursive MF is shown by (2.15), where represents the time series and the filter window has samples.

The recursive form of the MF presents an expressive noise attenuation capacity comparing with the nonrecursive form, this property leads the signal to a fast convergence [12].

The appropriate size of the MF can be defined comparing the variance of the signal filtered with the variance of the same signal filtered by an average filter. In this work it was chosen the exponentially weighted moving average (EWMA) filter [13], setting the factor properly. The EWMA filter is defined as follows: where is the estimated mean of the random variable at time and is the EWMA factor.

The absolute value of the difference between variances is defined as where is the size MF output and is the EWMA filter output related to factor.

2.4. -CUSUM Algorithm

The cumulative sum (CUSUM) algorithm is widely used to detect changes in the mean value of an independent Gaussian sequence. This algorithm is based on log-likelihood ratio and defined as follows [14]: where is the sufficient statistics, is probability density function with conditional parameter .

Before a change, the parameter is constant and equals to ; after the variation it assumes the value .

Considering (2.18) and a sampling set of size , the following decision function can be established: where the decision is stated as, and the hypothesis test is defined as, where is a suitable threshold.

The alarm time () after a change is defined by the following stopping rule: In this form, the CUSUM algorithm requires prior knowledge about the parameter and cannot be processed recursively. So, for online applications and considering the parameter as unknown, it is required a modified version of the CUSUM algorithm. It is presented in this work a modified version of the CUSUM based upon the weighted likelihood ratio and the sequential probability ratio test (SPRT) to suit the unknown and allow for the online processing as in [14]: where the stopping time (or alarm time) is defined as and is a weighting function.

Considering the prior knowledge about and of a Gaussian sequence with distribution concentrated on two points , and the weighted likelihood ratio becomes where is the signal-to-noise ratio (SNR), , and Solving (2.25), where Equation (2.27) is the probability ratio to noncentral parameter test of a distribution with one degree-of-freedom between values 0 and . So, by using this approach, we reach -CUSUM algorithm.

As a consequence of the previous derivation, is now defined as where By making a slight modification and considering the CUSUM algorithm as a repeated SPRT [14] with lower and upper thresholds fixed in 0 and , respectively, the recursive form of the -CUSUM algorithm takes the following form: where .

3. Fault Detection Algorithm

3.1. FD Algorithm Structure

The detection of faults in a IMU with redundant sensors can be performed by PV analysis. Under normal conditions, that is, bias and faults with null values, the PV should present a normal distribution . However, in actual systems, impulsive noise may occur. This situation is described in [7]. As stated above, the MF is used in order to remove those impulses and to improve the detection performed by -CUSUM algorithm. In Figure 2, it is shown the FD processing structure. Again, in this form, the algorithm only alarms the fault occurrence because it is not possible to define which sensor is under fault condition. After processing the signal from gyros by a data acquisition system, the data are filtered by a recursive MF (2.15) of size 3 or 5 in order to remove the spikes (or outliers) generated in this phase. In the sequence, the information (in /deg/s) is converted by a scale factor, SF (in deg/s). The SF block makes a polynomial conversion, whose degree is 7. The PV block performs the multiplication of the null space of the actual matrix (3.1) with the gyro output vector in deg/s, forming the parity vector. This PV is filtered, again, by another recursive MF of size 11. This size was obtained using (2.16) and (2.17), setting . The results obtained by processing (2.17) are shown in Figure 3. It can be seen that difference from the results obtained for sizes 9 to 15 is negligible. So, after some tests, it was concluded that size 11 is the best choice. The actual parameters used in the algorithm were obtained from the calibration of the IMU and are presented in [7]. The values and operations associated to blocks in the Figure 2 are given as follows.

For the actual sensor matrix (this matrix is obtained from calibration process) () and respective null space matrix () and generalized inverse (), For the scale factor block, the polynomials of degree 7 are, and the PV block processes the following operation: where is filtered by MF and applied to -CUSUM algorithm (2.31) in the form of .

The filtering block in the Figure 2 is an optional filter to comply with the GN&C block requirements.

3.2. Calibration of the FD Algorithm

The calibration of the FD algorithm was performed offline by using a time series extracted from the IMU movement on a 2DOF rotating table. The first parameter to be determined is the threshold (), obtained by use of the Kullback information [14] as follows: where is the average number of samples until the alarm time ().

Of course, the parameters and in (3.6) are unknown. So, to overcome this condition, it is assumed the SNR () equals to 1 and with properly size. In this work, it was considered equals to 60. Thus, the threshold is fixed as Once the threshold value is set and considered the actual values associated to the PV ( and ) as given in Table 1, it was performed the following computation the -CUSUM algorithm: The ratio in (3.8) is used to determine the actual value of that reach 0% of samples crossing the threshold in normal condition (without fault). In other words, 0% of false alarms. It is shown in Figure 4 the ratios for varying from 0.1 to 5, in whose range from 3.5 for SNR meets the objective. The efficiency of the calibration can be seen in Figures 5 and 6, it is shown the algorithm outputs for (SNR) and , respectively. The main information is presented in Figures 5 and 6(d), where it is illustrated the behavior of the decision function with respect to the threshold ().

4. Results

After the calibration phase, it was injected a step bias fault of magnitude 0.15 deg/s into one of sensors for the sample 10000. This fault generated a step variation in the PV according to Figures 7 and 8(a) and slopes as a result of the -CUSUM algorithm processing in the Figures 8(b)8(d). In the Figure 7(a), it is shown the efficiency of the median filter in impulsive noise reduction process, and in letter (b) is shown a comparison between MF and FExp in terms of filtering delay. In the decision function curve Figure 8(d), it is illustrated the difference between actual time occurrence of the fault () and the actual time alarm (), whose delay is 54 samples. In addition, there are another delays that should be computed which are associated to median filters (MF and MF). The total processing delay of the filters and algorithm is summarized in Table 2. In this table, it is presented the processing of the step faults with values 0.10, 0.15, 0.20, and 0.30 deg/s in order to compare the effect of the fault level on the detection delay. Considering the case of this work, where the FOGs are sampled at 100 Hz, the delay is less than 0.15 seconds for the fault level higher than 0.20 deg/s.

5. Conclusions

In this paper, a method based on -CUSUM algorithm combined with median filter was applied to detect faults in inertial measurement units with minimal redundancy of fiber optics gyros. The effectiveness of algorithm in case of low level step faults was demonstrated achieving the requirements of short delay to alarm with high reliability. In addition, the calibration technique presented here also proved to be efficient.