# IMU Fault Detection Based onχ2-CUSUM

Authors: Élcio Jeronimo de Oliveira; Hélio Koiti Kuga; Ijar Milagre da Fonseca
Journal: Mathematical Problems in Engineering
 (2012)

---

## 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 (χ2-CUSUM) only for fault detection (FD) applied to an IMU composed by four FOGs.


---

## Body

## 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 [1–4] 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 χ2-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 Figure1. 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 α=54.736°, the DCM is given by(2.1)H=(130631322-6613-22-66100).
The matrix H relates the sensor measurements (gi) with the angular rate in the main axes in the form of(2.2)go=Hω.
The estimate of the angular rate components in the main axes can be obtained from (2.2) in the following manner:(2.3)ω̂=H*go,(2.4)H*=(HTH)-1HT,
where go=[g1g2g3g4]T is the vector of the gyro outputs, ω=[ωxωyωz]T is the angular rate vector about main axes, ω̂ is the angular rate estimate vector, and H* is the generalized inverse of H.Tetrahedral base.
(a)(b)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:(2.5)g̃o=Hω+δgo+f+ηs,
where δgo is the bias vector whose magnitude is constant, f is the fault vector, and ηs is a random term vector (Gaussian noise). Applying the singular value decomposition (SVD) on H, 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 H as follows:(2.6)UTHV=Λ=(𝚺0),(2.7)H=UΛVT,
where U, Λ, and VT are matrices obtained from SVD of H. The matrix Σ is a diagonal matrix whose elements are the eigenvalues of H. The superscript (T) indicates the transpose. Applying (2.7) into (2.5) and multiplying both sides by UT, it can be obtained the following relationship:(2.8)UTg̃o=ΛVTω+UT(δgo+f+ηs).
Partitioning U as follows:(2.9)U=[U1⋮U2],
where U1∈ℝ4×3 and U2∈ℝ4×1, and applying it into (2.8), the resulting equations are(2.10)U1Tg̃o=(𝚺V)ω+U1T(δgo+f+ηs),(2.11)U2Tg̃o=U2T(δgo+f+ηs)=p,(2.12)p=Cg̃o.
Equation (2.10) leads to least squares estimate of ω, what is equivalent to (2.3), and can be expressed by(2.13)ω̂LS=(𝚺V)-1U1Tg̃o.
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, U2T (or C) is null space of H and p 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 (k), shift the window to process the time (k+1) 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 (k), replace the original sample x(k) by the filter output y(k), 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),(2.14)y(k)=median{x(k-N),…,x(k),…,x(k+N)},(2.15)y(k)=median{y(k-N),…,y(k-1),x(k),…,x(k+N)},
where x represents the time series and the filter window has 2N+1 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: (2.16)μ̅n=βμ̅n-1+(1-β)xn,
where μ̅ is the estimated mean of the random variable x at time n and β is the EWMA factor.The absolute value of the difference between variances is defined as(2.17)Δvar(j)=abs[var(MF(j))-var(FExp(β))],
where MF(j) is the size j MF output and FExp(β) is the EWMA filter output related to β factor.


### 2.4.χ2-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]:(2.18)s(y)=ln⁡pθ1(y)pθ0(y),
where s(y) is the sufficient statistics, pθ(y) is probability density function with conditional parameter θ.Before a change, the parameterθ is constant and equals to θ0; after the variation it assumes the value θ1.Considering (2.18) and a sampling set of size N, the following decision function can be established:(2.19)Sjk=∑i=jksi,si=ln⁡pθ1(yi)pθ0(yi),
where the decision is stated as,(2.20)d={0ifS1N<λ,H0ischosen1ifS1N≥λ,H1ischosen
and the hypothesis test is defined as,(2.21)H0:θ=θ0,H1:θ=θ1,
where λ is a suitable threshold.The alarm time (ta) after a change is defined by the following stopping rule:(2.22)ta=N⋅min⁡{K:dK=1}.
In this form, the CUSUM algorithm requires prior knowledge about the parameter θ1 and cannot be processed recursively. So, for online applications and considering the parameter θ1 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 θ1 and allow for the online processing as in [14]: (2.23)Λ̃jk=∫-∞∞pθ1(yj,…,yk)pθ0(yj,…,yk)dF(θ1),
where the stopping time (or alarm time) ta is defined as(2.24)ta=min⁡⁡{k:max⁡1≤j≤kln⁡⁡Λ̃jk≥λ}
and dF(θ1) is a weighting function.Considering the prior knowledge aboutσ2 and μ0 of a Gaussian sequence with distribution F(θ1)=F(μ) concentrated on two points μ0-ν, and μ0+ν the weighted likelihood ratio becomes(2.25)Λ̃jk=∫-∞∞exp⁡[bS̃jk-b22(k-j+1)]dF(ν),
where b=ν/σ is the signal-to-noise ratio (SNR), ν=μ1-μ0, and(2.26)S̃jk=1σ∑i=jk(yi-μ0).
Solving (2.25),(2.27)Λ̃jk=cosh⁡(bS̃jk)e-(b2/2)(k-j+1)=cosh⁡[b(k-j+1)χjk]e-(b2/2)(k-j+1),
where(2.28)χjk=1k-j+1|S̃jk|.
Equation (2.27) is the probability ratio to noncentral parameter test of a χ2 distribution with one degree-of-freedom between values 0 and (k-j+1)b2. So, by using this approach, we reach χ2-CUSUM algorithm.As a consequence of the previous derivation,ta is now defined as(2.29)ta=min⁡{k:gk≥λ},
where(2.30)gk=max⁡1≤j≤k[ln⁡⁡cosh⁡⁡(bS̃jk)-b22(k-j+1)].
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 χ2-CUSUM algorithm takes the following form:(2.31)gk=(S̆k-Nk+1k)+,S̆k-Nk+1k=-12Nkb2+ln⁡⁡cosh⁡⁡(bS̃k-Nk+1k),S̃k-Nk+1k=S̅k,S̅k=S̅k-11{gk-1>0}+yk-μ0σ,
where Nk=Nk-11{gk-1>0}+1.


## 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 Figure1. 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 α=54.736°, the DCM is given by(2.1)H=(130631322-6613-22-66100).
The matrix H relates the sensor measurements (gi) with the angular rate in the main axes in the form of(2.2)go=Hω.
The estimate of the angular rate components in the main axes can be obtained from (2.2) in the following manner:(2.3)ω̂=H*go,(2.4)H*=(HTH)-1HT,
where go=[g1g2g3g4]T is the vector of the gyro outputs, ω=[ωxωyωz]T is the angular rate vector about main axes, ω̂ is the angular rate estimate vector, and H* is the generalized inverse of H.Tetrahedral base.
(a)(b)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:(2.5)g̃o=Hω+δgo+f+ηs,
where δgo is the bias vector whose magnitude is constant, f is the fault vector, and ηs is a random term vector (Gaussian noise). Applying the singular value decomposition (SVD) on H, 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 H as follows:(2.6)UTHV=Λ=(𝚺0),(2.7)H=UΛVT,
where U, Λ, and VT are matrices obtained from SVD of H. The matrix Σ is a diagonal matrix whose elements are the eigenvalues of H. The superscript (T) indicates the transpose. Applying (2.7) into (2.5) and multiplying both sides by UT, it can be obtained the following relationship:(2.8)UTg̃o=ΛVTω+UT(δgo+f+ηs).
Partitioning U as follows:(2.9)U=[U1⋮U2],
where U1∈ℝ4×3 and U2∈ℝ4×1, and applying it into (2.8), the resulting equations are(2.10)U1Tg̃o=(𝚺V)ω+U1T(δgo+f+ηs),(2.11)U2Tg̃o=U2T(δgo+f+ηs)=p,(2.12)p=Cg̃o.
Equation (2.10) leads to least squares estimate of ω, what is equivalent to (2.3), and can be expressed by(2.13)ω̂LS=(𝚺V)-1U1Tg̃o.
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, U2T (or C) is null space of H and p 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 (k), shift the window to process the time (k+1) 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 (k), replace the original sample x(k) by the filter output y(k), 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),(2.14)y(k)=median{x(k-N),…,x(k),…,x(k+N)},(2.15)y(k)=median{y(k-N),…,y(k-1),x(k),…,x(k+N)},
where x represents the time series and the filter window has 2N+1 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: (2.16)μ̅n=βμ̅n-1+(1-β)xn,
where μ̅ is the estimated mean of the random variable x at time n and β is the EWMA factor.The absolute value of the difference between variances is defined as(2.17)Δvar(j)=abs[var(MF(j))-var(FExp(β))],
where MF(j) is the size j MF output and FExp(β) is the EWMA filter output related to β factor.


## 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 (k), shift the window to process the time (k+1) 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 (k), replace the original sample x(k) by the filter output y(k), 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),(2.14)y(k)=median{x(k-N),…,x(k),…,x(k+N)},(2.15)y(k)=median{y(k-N),…,y(k-1),x(k),…,x(k+N)},
where x represents the time series and the filter window has 2N+1 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: (2.16)μ̅n=βμ̅n-1+(1-β)xn,
where μ̅ is the estimated mean of the random variable x at time n and β is the EWMA factor.The absolute value of the difference between variances is defined as(2.17)Δvar(j)=abs[var(MF(j))-var(FExp(β))],
where MF(j) is the size j MF output and FExp(β) is the EWMA filter output related to β factor.


## 2.4.χ2-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]:(2.18)s(y)=ln⁡pθ1(y)pθ0(y),
where s(y) is the sufficient statistics, pθ(y) is probability density function with conditional parameter θ.Before a change, the parameterθ is constant and equals to θ0; after the variation it assumes the value θ1.Considering (2.18) and a sampling set of size N, the following decision function can be established:(2.19)Sjk=∑i=jksi,si=ln⁡pθ1(yi)pθ0(yi),
where the decision is stated as,(2.20)d={0ifS1N<λ,H0ischosen1ifS1N≥λ,H1ischosen
and the hypothesis test is defined as,(2.21)H0:θ=θ0,H1:θ=θ1,
where λ is a suitable threshold.The alarm time (ta) after a change is defined by the following stopping rule:(2.22)ta=N⋅min⁡{K:dK=1}.
In this form, the CUSUM algorithm requires prior knowledge about the parameter θ1 and cannot be processed recursively. So, for online applications and considering the parameter θ1 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 θ1 and allow for the online processing as in [14]: (2.23)Λ̃jk=∫-∞∞pθ1(yj,…,yk)pθ0(yj,…,yk)dF(θ1),
where the stopping time (or alarm time) ta is defined as(2.24)ta=min⁡⁡{k:max⁡1≤j≤kln⁡⁡Λ̃jk≥λ}
and dF(θ1) is a weighting function.Considering the prior knowledge aboutσ2 and μ0 of a Gaussian sequence with distribution F(θ1)=F(μ) concentrated on two points μ0-ν, and μ0+ν the weighted likelihood ratio becomes(2.25)Λ̃jk=∫-∞∞exp⁡[bS̃jk-b22(k-j+1)]dF(ν),
where b=ν/σ is the signal-to-noise ratio (SNR), ν=μ1-μ0, and(2.26)S̃jk=1σ∑i=jk(yi-μ0).
Solving (2.25),(2.27)Λ̃jk=cosh⁡(bS̃jk)e-(b2/2)(k-j+1)=cosh⁡[b(k-j+1)χjk]e-(b2/2)(k-j+1),
where(2.28)χjk=1k-j+1|S̃jk|.
Equation (2.27) is the probability ratio to noncentral parameter test of a χ2 distribution with one degree-of-freedom between values 0 and (k-j+1)b2. So, by using this approach, we reach χ2-CUSUM algorithm.As a consequence of the previous derivation,ta is now defined as(2.29)ta=min⁡{k:gk≥λ},
where(2.30)gk=max⁡1≤j≤k[ln⁡⁡cosh⁡⁡(bS̃jk)-b22(k-j+1)].
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 χ2-CUSUM algorithm takes the following form:(2.31)gk=(S̆k-Nk+1k)+,S̆k-Nk+1k=-12Nkb2+ln⁡⁡cosh⁡⁡(bS̃k-Nk+1k),S̃k-Nk+1k=S̅k,S̅k=S̅k-11{gk-1>0}+yk-μ0σ,
where Nk=Nk-11{gk-1>0}+1.


## 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𝒩(0,σ2). 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 χ2-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 mV/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 H̃ 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 β=0.95. 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.Figure 2
Gyro processing and fault detection block diagram. MF: median filter (MF(1): size 3//MF(2): size 11); SF: scale factor; GN&C: guidance, navigation, and control; PV: parity vector.Figure 3
Absolute value of the difference between variances of theFM(j) and FExp(0.95).For the actual sensor matrix (this matrix is obtained from calibration process) (H̃) and respective null space matrix (C̃) and generalized inverse (H*), (3.1)H̃=(0.57868624-0.004324360.815538790.577337660.70734792-0.407848170.57695344-0.70679030-0.409355820.999998050.001009950.00169158),(3.2)C̃=[-0.40955-0.40903-0.405850.70729],C̃H̃=0,H*=(H̃TH̃)-1H̃T.
For the scale factor block, the polynomials of degree 7 are,(3.3)go1(°s)=0.0657gv7-3.64×10-4gv6-7.799×10-2gv5+9.911×10-4gv4+1.200gv3+3.085×10-4gv2+37.432gv-3.724×10-3,go2(°s)=0.09452gv7+1.083×10-5gv6-0.1965gv5-2.041×10-3gv4+1.5219gv3+6.910×10-3gv2+38.305gv-4.615×10-3,go3(°s)=0.0761gv7-2.898×10-4gv6-0.13277gv5-8.253×10-4gv4+1.3711gv3+2.326×10-3gv2+38.234gv-1.1713×10-2,go4(°s)=0.0814gv7-1.0544×10-4gv6-0.18366gv5-1.3751×10-3gv4+1.4068gv3+3.9596×10-3gv2+37.282gv-4.1806×10-3,
and the PV block processes the following operation:(3.4)p=C̃g̃o,
where p is filtered by MF(2) and applied to χ2-CUSUM algorithm (2.31) in the form of yk.The filtering block in the Figure2 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:(3.5)l̅≈λK(ν),K(ν)=ν22σ2=12b2,(3.6)λ≈12b2l̅,
where l̅ is the average number of samples until the alarm time (ta).Of course, the parametersb and l̅ in (3.6) are unknown. So, to overcome this condition, it is assumed the SNR (b) equals to 1 and l̅ with properly size. In this work, it was considered l̅ equals to 60. Thus, the threshold is fixed as(3.7)λ=1212×60=30.
Once the threshold value is set and considered the actual values associated to the PV (σ and μ0) as given in Table 1, it was performed the following computation the χ2-CUSUM algorithm:(3.8)%=100×Totalofgk>λTotalnumberofsamples.
The ratio in (3.8) is used to determine the actual value of b 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 b 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) b=0.1 and b=3.5, 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 (λ=30).Table 1
Tuning parameters forχ2-CUSUM algorithm.
σ(°/s)μ0(°/s)bλ0.0220.0020.1 to 530Figure 4
Ratios forb varying from 0.1 to 5.χ2-CUSUM results for 20000 samples of a movement sequence without faults with (SNR) b=0.1. (a) PV; (b) S̅k; (c) S̆k; (d) gk.
(a)(b)(c)(d)χ2-CUSUM results for 20000 samples of a movement sequence without faults with (SNR) b=3.5. (a) PV; (b) S̅k; (c) S̆k; (d) gk.
(a)(b)(c)(d)


## 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𝒩(0,σ2). 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 χ2-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 mV/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 H̃ 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 β=0.95. 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.Figure 2
Gyro processing and fault detection block diagram. MF: median filter (MF(1): size 3//MF(2): size 11); SF: scale factor; GN&C: guidance, navigation, and control; PV: parity vector.Figure 3
Absolute value of the difference between variances of theFM(j) and FExp(0.95).For the actual sensor matrix (this matrix is obtained from calibration process) (H̃) and respective null space matrix (C̃) and generalized inverse (H*), (3.1)H̃=(0.57868624-0.004324360.815538790.577337660.70734792-0.407848170.57695344-0.70679030-0.409355820.999998050.001009950.00169158),(3.2)C̃=[-0.40955-0.40903-0.405850.70729],C̃H̃=0,H*=(H̃TH̃)-1H̃T.
For the scale factor block, the polynomials of degree 7 are,(3.3)go1(°s)=0.0657gv7-3.64×10-4gv6-7.799×10-2gv5+9.911×10-4gv4+1.200gv3+3.085×10-4gv2+37.432gv-3.724×10-3,go2(°s)=0.09452gv7+1.083×10-5gv6-0.1965gv5-2.041×10-3gv4+1.5219gv3+6.910×10-3gv2+38.305gv-4.615×10-3,go3(°s)=0.0761gv7-2.898×10-4gv6-0.13277gv5-8.253×10-4gv4+1.3711gv3+2.326×10-3gv2+38.234gv-1.1713×10-2,go4(°s)=0.0814gv7-1.0544×10-4gv6-0.18366gv5-1.3751×10-3gv4+1.4068gv3+3.9596×10-3gv2+37.282gv-4.1806×10-3,
and the PV block processes the following operation:(3.4)p=C̃g̃o,
where p is filtered by MF(2) and applied to χ2-CUSUM algorithm (2.31) in the form of yk.The filtering block in the Figure2 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:(3.5)l̅≈λK(ν),K(ν)=ν22σ2=12b2,(3.6)λ≈12b2l̅,
where l̅ is the average number of samples until the alarm time (ta).Of course, the parametersb and l̅ in (3.6) are unknown. So, to overcome this condition, it is assumed the SNR (b) equals to 1 and l̅ with properly size. In this work, it was considered l̅ equals to 60. Thus, the threshold is fixed as(3.7)λ=1212×60=30.
Once the threshold value is set and considered the actual values associated to the PV (σ and μ0) as given in Table 1, it was performed the following computation the χ2-CUSUM algorithm:(3.8)%=100×Totalofgk>λTotalnumberofsamples.
The ratio in (3.8) is used to determine the actual value of b 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 b 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) b=0.1 and b=3.5, 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 (λ=30).Table 1
Tuning parameters forχ2-CUSUM algorithm.
σ(°/s)μ0(°/s)bλ0.0220.0020.1 to 530Figure 4
Ratios forb varying from 0.1 to 5.χ2-CUSUM results for 20000 samples of a movement sequence without faults with (SNR) b=0.1. (a) PV; (b) S̅k; (c) S̆k; (d) gk.
(a)(b)(c)(d)χ2-CUSUM results for 20000 samples of a movement sequence without faults with (SNR) b=3.5. (a) PV; (b) S̅k; (c) S̆k; (d) gk.
(a)(b)(c)(d)


## 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 Figures7 and 8(a) and slopes as a result of the χ2-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 (t0=10000) and the actual time alarm (ta=10054), whose delay is 54 samples. In addition, there are another delays that should be computed which are associated to median filters (MF(1) and MF(2)). 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.Table 2
Fault detection delay as a function of fault magnitude given in number of samples with respect to the parametersb=3.5 and λ=30.
Step fault1st MF delay2nd MF delayDetectionTotal delay(°/s)(MF3)(MF11)delay(ta-t0)0.10154304360.151554600.20159150.3015410PV with step fault for sample 10000. (a) PV nonfiltered (red) and filtered by MF of size 11 (blue); (b) PV filtered by MF of size 11 (blue) and by FExp withβ=0.95 (red); the zoom in (b) indicates the difference between MF and FExp during the transition time.
(a)(b)χ2-CUSUM results for 20000 samples of a movement sequence with a step fault at sample 10000 and magnitude 0.15°/s; (SNR) b=3.5. (a) PV; (b) S̅k; (c) S̆k; (d) gk.
(a)(b)(c)(d)


## 5. Conclusions

In this paper, a method based onχ2-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.



---

*Source: 740752-2012-04-11.xml*