# Fault Detection and Isolation in Inertial Measurement Units Based onχ2-CUSUM and Wavelet Packet

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

---

## Abstract

The aim of this paper is to present a fault detection algorithm (FDI) based on signal
processing techniques developed for an inertial measurement unit (IMU) with minimal
redundancy of fiber optic gyros. In this work   the recursive median filter is applied in
order to remove impulses (outliers) arising from data acquisition process and parity vector
operations, improving the fault detection and isolation performance. The FDI algorithm is
divided into two blocks: fault detection (FD) and fault isolation (FI). The FD part of the
algorithm is used to guarantee the reliability of the isolation part and is based on parity
vector analysis usingχ2-CUSUM algorithm. The FI part is performed using parity space
projection of the energy subbands obtained from wavelet packet decomposition. This
projection is an extension of clustering analysis based on singular value decomposition
(SVD) and principal component analysis (PCA). The results of the FD and FI algorithms
have shown the effectiveness of the proposed method, in which the FD algorithm is capable
of indicating  the low-level  step bias fault with short delay and a high index of correct decisions
of the FI algorithm also with low-level step bias fault.


---

## Body

## 1. Introduction

The main goal of fault detection, isolation, and recovery (FDIR) is to effectively detect faults and accurately isolate them to a failed component in the shortest time possible. This capability leads to reduction in diagnostic time or downtime in general, increasing the system availability. A good inherent diagnostic of a system also enhances the self-confidence in operating system, the main aspect of mission success. FDIR is especially required to an on-orbit system where maintenance may be difficult, very expensive, or even impossible. The spacecraft AOCS (Attitude and Orbit Control Subsystem) includes components such as sensors and actuators. Among the attitude control devices, the gyros are widely used for the sake of attitude control, mainly when high-accuracy requirements are imposed on the AOCS. This work belongs to this area and is addressed in the scenario of the orbital dynamics for spacecraft AOCS including inertial measurement unities (IMU) and FDIR software. This paper does not cover the recovery problem associated with FDIR but only with the FDI, that is, the fault detection and isolation problem. The methodology is based on signal processing techniques developed for an inertial measurement unit (IMU) with minimal redundancy of fiber optic gyros.In the FDI problem with minimal redundancy, we are limited only to detect faults, generally with the aid of the parity vector analysis [1–6]. For minimal redundancy, we mean the use of four inertial sensors (gyros or accelerometers) to measure a tridimensional space where, at least, one of the sensors should be skewed [7]. The identification of the faulty sensor in the minimal redundancy case is usually performed by signal processing techniques [8]. In the last years, the increase of the wavelets applications for fault (or variation) detection [8–11] and classification [12] has demonstrated the potential of the wavelet for IMU status monitoring. In addition, actual digital systems, mainly in prototype stage, generate high quantity of impulsive noise (outliers) as a consequence of the data acquisition system and/or parity vector operations [13]. To overcome this problem, removing the outliers without degenerating important variations in the sensor signals, the median filter (MF) is used in the recursive form. This kind of filter, widely used in image processing, is applied here because it preserves the edge transitions without distortions [14] and is a robust estimator [15]. These properties are helpful in the wavelet analyses to preserve important information in the details. Another advantage in the MF use is the increase of the performance in the change detection [13]. In the pattern recognition field, the singular value decomposition (SVD) and the principal component analysis (PCA) are useful tools in the practical approach to detect hidden structure in the form of clustered projections in noisy data and with reduced number of variables [16]. Therefore, this paper is dedicated to develop an algorithm applied to faulty sensor isolation (FI) in IMU constructed with four fiber optic gyros (FOG) in a complementary development of the fault detection algorithm (FD) discussed in detail by de Oliveira [6]. The FD algorithm is based on parity vector analysis processed by χ2-CUSUM algorithm [17] for change detections and is designed in order to guarantee the reliability of the isolation part. The FI part is designed using an approach based on parity space projection of the energy subbands obtained from wavelet packet decomposition. This projection is an extension of clustering analysis based on singular value decomposition (SVD) and principal component analysis (PCA).This paper is organized as follows. Section2 discusses the theoretical bases applied in this work composed by geometrical and parity vector (PV) analysis, recursive median filter (MF), χ2-CUSUM algorithm, wavelet packet (WP) decomposition, singular value decomposition (SVD), and principal component analysis (PCA). Section 3 shows the development of the FDI algorithm blocks and their relationships. In that section, the FI model is derived with the parity space projection of the energy subbands obtained by WP decomposition. In Section 4, the results and respective parameters used for FI algorithms are discussed. Section 5 presents an  overview of the work and conclusions. The appendix  shows the parameters obtained by FI algorithm calibration.


## 2. Background

The theoretical background and calibration methodology applied for FD algorithm (geometrical and parity vector analyses, recursive median filter, andχ2-CUSUM algorithm) were completely developed in previous work [6], and here they will be omitted. This work is dedicated to design a complete FDI algorithm, and, to reach this goal, it is needed to introduce a new set of mathematical tools to develop the FI algorithm. Therefore, in this section, the theoretical aspects related to geometry, WP, SVD, and PCA will be developed.


### 2.1. Geometry

The geometrical aspect plays an important role in the inertial measurement units design. Several works [1, 3, 7, 18] address the matter of geometrical analyses, and mainly in the optimal angular relationship between sensors. In this work, a minimal redundancy (four gyros) configuration of fiber optic gyros is used as shown in Figure 1. The nominal sensor matrix (direct cosine matrix) related to this configuration is expressed by (1) and measurement equation by (2):
(1)H=(130631322-6613-22-66100)(2)g=Hω+ν+b+f,Tetrahedral base and respective gyro axes orientation.
(a)
(b)whereg is the sensor measurement vector (m×1); H is the observation matrix (m×3) relating the sensor axes with a triorthogonal axes system; ω is the state vector (3×1) in the tri-orthogonal axes system (ω=[ωxωyωz]T); ν is the Gaussian white noise vector (m×1); and  b is the bias vector (m×1); f is the fault vector (m×1); m is the number of sensors.The parity vector (p) associated with this configuration is given by [13]
(3)Cg=C(Hω+ν+b+f),p=C(Hω)+C(ν+b+f),where the matrixC is obtained from the null space of H, so that
(4)CH=0.In the absence of biases and faults, the parity vector is a white Gaussian noise vector, and in this configuration (four sensors in a skewed arrangement), it is possible to detect only the fault occurrence.


### 2.2. Wavelet Packet

The wavelet packet (WP) decomposition is a method based on conventional wavelet transform theory whose difference is the processing of both filter outputs. In this method, the low-pass and high-pass filter outputs are filtered at each level of decomposition in the form of a binary symmetric decomposition tree. This method is useful when a high spectral resolution analysis with low computational cost is needed. The WP decomposition structure is shown in Figure2. At each level of decomposition, the WP generates 2j blocks with dim(s)/2j coefficients, where j is the decomposition level and dim(s) is the size of the buffer (number of samples to be processed). The downsampling operator (↓2) follows the LP and HP filters at each level in order to eliminate the redundant information [19].Figure 2
Wavelet packet decomposition tree up to levelj=3 from buffer with 256 samples. yi: signal samples; yi(): coefficients after downsampling at each level of decomposition; LP: low-pass filter output (a-approximations); HP: high-pass filter output (d-details).


### 2.3. SVD and PCA

The singular value decomposition (SVD) is a mathematical operation defined in the linear algebra branch, that allows the real or complex matrix factorization. Its application is very useful in signal processing and analyses [4, 16, 20]. The SVD is defined according, Theorem 1 [21].Theorem 1 (SVD).
LetA be a matrix (real or complex) m×n of rank r. Therefore, there exist an unitary matrix U (orthogonal if A is real) with order m, an unitary matrix V (orthogonal) with order n, and a diagonal matrix S with order r (with positive values), so that,
(5)A=U[S000]V*
or(6)A=UΣV*,
where,V* is the conjugate transpose of  V and Σ is m×n. The values in the main diagonal of  S are the singular values of A.The SVD and the eigenvalues and eigenvectors decomposition are related as follows.(i)
The left-hand side singular vectors ofA (columns of U) are the eigenvectors of AA*.(ii)
The right-hand side singular vectors ofA (columns of V) are the eigenvectors of A*A.(iii)
The nonnull singular values ofS are the square root of the non-null eigenvectors of AA* or A*A.Being the covariance matrix ofA given by A*A, the singular values of A are the square root of the eigenvalues of its covariance matrix.
(7)A*A=VS2V*.By using SVD the principal components from a given matrix can be obtained. Mathematically, the PCA is defined as a linear operation that projects a data set into a new coordinate systems, so that the major variance component be associated with the first coordinate, the second major variance component with the second coordinate, and so on [22]. In this way, a set of correlated variables are transformed into another uncorrelated variable set called principal components. The principal components can be obtained by decomposition of eigenvalues from covariance matrix or using SVD. The main advantage of use the PCA is the possibility of reducing the dimensionality of the data set space (size n), by selecting the m main components from the space of size n (m<n), ensuring an optimal orthogonal transformation [22]. Mathematically, PCA can be obtained with the aid of the SVD by defining a set of data (X) arranged as a matrix according to Table 1, where they are stacked J sequences of K samples, being r the rank of the matrix. In order to calculate the effective direction of the maximum variation, it is extracted the average value from each line as follows:
(8)x-j=1K∑k=1K‍xj,k,j=1⋯J,MX=X-x-1K,1K=[11⋯1]1×K.Table 1
Block of data set with rankn.
Seq.
Samples
1
x
11
x
12
⋯
x
1
K
2
x
21
x
22
⋯
x
2
K
3
x
31
x
32
⋯
x
3
K
⋮
⋮
J
x
J
1
x
J
2
⋯
x
J
KWith data set centered at zero mean, the SVD on the matrixMX can be performed  as follows:
(9)[USV]=svd(MX).By sorting the data in decreasing order, the major values, and respective eigenvectors, can be selected considering a properly threshold criteria. By definingES as a pseudoenergy related to sum of all singular values obtained in the form,
(10)ES=∑j=1J‍sjj,wheresjj are the diagonal elements of S; and the cumulative pseudoenergy es(j) as a cumulative the sum of the pseudoenergies related to each eigenvector:
(11)es(j)=∑p=1j‍spp,j=1,…,J,a threshold is established in order to define a numberL of principal components (eigenvectors) that should be selected as follows:
(12)es(j=L)ES≥th.After defining theL most significant components, the matrix V can be resized considering only the first L columns. Then,
(13)W=Vpq,wherep=1,…,J and q=1,…,L.


## 2.1. Geometry

The geometrical aspect plays an important role in the inertial measurement units design. Several works [1, 3, 7, 18] address the matter of geometrical analyses, and mainly in the optimal angular relationship between sensors. In this work, a minimal redundancy (four gyros) configuration of fiber optic gyros is used as shown in Figure 1. The nominal sensor matrix (direct cosine matrix) related to this configuration is expressed by (1) and measurement equation by (2):
(1)H=(130631322-6613-22-66100)(2)g=Hω+ν+b+f,Tetrahedral base and respective gyro axes orientation.
(a)
(b)whereg is the sensor measurement vector (m×1); H is the observation matrix (m×3) relating the sensor axes with a triorthogonal axes system; ω is the state vector (3×1) in the tri-orthogonal axes system (ω=[ωxωyωz]T); ν is the Gaussian white noise vector (m×1); and  b is the bias vector (m×1); f is the fault vector (m×1); m is the number of sensors.The parity vector (p) associated with this configuration is given by [13]
(3)Cg=C(Hω+ν+b+f),p=C(Hω)+C(ν+b+f),where the matrixC is obtained from the null space of H, so that
(4)CH=0.In the absence of biases and faults, the parity vector is a white Gaussian noise vector, and in this configuration (four sensors in a skewed arrangement), it is possible to detect only the fault occurrence.


## 2.2. Wavelet Packet

The wavelet packet (WP) decomposition is a method based on conventional wavelet transform theory whose difference is the processing of both filter outputs. In this method, the low-pass and high-pass filter outputs are filtered at each level of decomposition in the form of a binary symmetric decomposition tree. This method is useful when a high spectral resolution analysis with low computational cost is needed. The WP decomposition structure is shown in Figure2. At each level of decomposition, the WP generates 2j blocks with dim(s)/2j coefficients, where j is the decomposition level and dim(s) is the size of the buffer (number of samples to be processed). The downsampling operator (↓2) follows the LP and HP filters at each level in order to eliminate the redundant information [19].Figure 2
Wavelet packet decomposition tree up to levelj=3 from buffer with 256 samples. yi: signal samples; yi(): coefficients after downsampling at each level of decomposition; LP: low-pass filter output (a-approximations); HP: high-pass filter output (d-details).


## 2.3. SVD and PCA

The singular value decomposition (SVD) is a mathematical operation defined in the linear algebra branch, that allows the real or complex matrix factorization. Its application is very useful in signal processing and analyses [4, 16, 20]. The SVD is defined according, Theorem 1 [21].Theorem 1 (SVD).
LetA be a matrix (real or complex) m×n of rank r. Therefore, there exist an unitary matrix U (orthogonal if A is real) with order m, an unitary matrix V (orthogonal) with order n, and a diagonal matrix S with order r (with positive values), so that,
(5)A=U[S000]V*
or(6)A=UΣV*,
where,V* is the conjugate transpose of  V and Σ is m×n. The values in the main diagonal of  S are the singular values of A.The SVD and the eigenvalues and eigenvectors decomposition are related as follows.(i)
The left-hand side singular vectors ofA (columns of U) are the eigenvectors of AA*.(ii)
The right-hand side singular vectors ofA (columns of V) are the eigenvectors of A*A.(iii)
The nonnull singular values ofS are the square root of the non-null eigenvectors of AA* or A*A.Being the covariance matrix ofA given by A*A, the singular values of A are the square root of the eigenvalues of its covariance matrix.
(7)A*A=VS2V*.By using SVD the principal components from a given matrix can be obtained. Mathematically, the PCA is defined as a linear operation that projects a data set into a new coordinate systems, so that the major variance component be associated with the first coordinate, the second major variance component with the second coordinate, and so on [22]. In this way, a set of correlated variables are transformed into another uncorrelated variable set called principal components. The principal components can be obtained by decomposition of eigenvalues from covariance matrix or using SVD. The main advantage of use the PCA is the possibility of reducing the dimensionality of the data set space (size n), by selecting the m main components from the space of size n (m<n), ensuring an optimal orthogonal transformation [22]. Mathematically, PCA can be obtained with the aid of the SVD by defining a set of data (X) arranged as a matrix according to Table 1, where they are stacked J sequences of K samples, being r the rank of the matrix. In order to calculate the effective direction of the maximum variation, it is extracted the average value from each line as follows:
(8)x-j=1K∑k=1K‍xj,k,j=1⋯J,MX=X-x-1K,1K=[11⋯1]1×K.Table 1
Block of data set with rankn.
Seq.
Samples
1
x
11
x
12
⋯
x
1
K
2
x
21
x
22
⋯
x
2
K
3
x
31
x
32
⋯
x
3
K
⋮
⋮
J
x
J
1
x
J
2
⋯
x
J
KWith data set centered at zero mean, the SVD on the matrixMX can be performed  as follows:
(9)[USV]=svd(MX).By sorting the data in decreasing order, the major values, and respective eigenvectors, can be selected considering a properly threshold criteria. By definingES as a pseudoenergy related to sum of all singular values obtained in the form,
(10)ES=∑j=1J‍sjj,wheresjj are the diagonal elements of S; and the cumulative pseudoenergy es(j) as a cumulative the sum of the pseudoenergies related to each eigenvector:
(11)es(j)=∑p=1j‍spp,j=1,…,J,a threshold is established in order to define a numberL of principal components (eigenvectors) that should be selected as follows:
(12)es(j=L)ES≥th.After defining theL most significant components, the matrix V can be resized considering only the first L columns. Then,
(13)W=Vpq,wherep=1,…,J and q=1,…,L.


## 3. FDI Algorithm

The basic structure of the gyros processing with FDI blocks embedded is shown in Figure3. In that structure, each gyro output is filtered by a recursive median filter (MF-1) in order to remove the outliers (impulsive noise), and in the sequence the gyro measurements are adjusted properly by scale factor (SF) block. After that, the measurements follow two different branches. The first one points towards GN&C (guidance, navigation, and control) block through a generic filter and sensor matrix blocks. As a consequence of the geometry presented in Figure 1, and respective model discussed through (1)-(2), the sensor matrix block performs the transformation from sensors space (tetrahedron) to state space (rotations about x, y, and z) according to (14), which represents a Least Square Estimation of the state, considering sensor measurements with errors. The second branch points to FDI blocks and it is divided in to two distinct sets. One of them performs the on-line fault detection processing, and the other performs the faulty sensor identification (fault isolation) after a fault is detected. Once the faulty sensor is identified, its influence is removed from the sensor matrix (H*), and the GN&C starts to operate only with three  pieces of FOG information:
(14)ω^=H*g(15)H*=(HTH)-1HT.Figure 3
Connection between Scale Factor block and FD/FI blocks is missed.


### 3.1. Fault Detection Algorithm

The detection of faults in an IMU with redundant sensors can be performed by PV analysis as defined in (3). 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 [13]. In the FD block, the recursive MF is used in order to remove impulsive noise and to improve the detection performed by χ2-CUSUM algorithm. In Figure 3 the FD processing structure in the left branch is shown. Under minimal redundancy configuration (four gyros), the FD algorithm can only indicate the fault occurrence. So, this information is used to start the identification process with high confidence of the fault occurrence. After the signal from gyros is processed by PV block according (3), the measurements are filtered by a second recursive MF (MF-2) to remove the spikes (or outliers) generated in the PV processing [6]. The filtered PV is applied to χ2-CUSUM block that performs the fault detection according to the algorithm and calibration presented in previous work [6].


### 3.2. Fault Isolation Algorithm

The fault isolation (FI) algorithm begins processing the data information stored in the buffers when the fault occurrence is indicated by FD block. The four buffers (one for each gyro) operate in FIFO mode. The operation of the FI algorithm is based on the analysis of the principal components in the energy subbands of the data stored in the buffers. The energy subbands are obtained by using the wavelet packet decomposition up to certain level. The idea behind this approach is to obtain asignature of the sensors in normal conditions given in terms of null space projection matrix and then use this matrix as kernel of FI analysis as demonstrated in the sequence.


#### 3.2.1. FI Algorithm Model

The modeling of the FI algorithm begins with WP decomposition of a large amount of data from sensors divided in to pieces with properly defined size. The WP decomposition follows the schema illustrated in Figure2. Once the appropriated level of decomposition is reached, the respective subband elements are processed according to (16) that represents the average pseudo-energy level at each subband:
(16)ey=1N(ysb)T(ysb),whereN and ysb are the amounts of wavelet coefficients and respective wavelet coefficients vector at given subband.At each decomposition levelj, we have 2j subbands. The first subband (superscripts a,aa,aaa,…) is composed by approximation coefficients and represents the sensor state information (low-pass filter output), and this information will not be used in FI processing. The remaining subbands will form the pseudo-energy matrix (j≥2) because those subbands contain the noise information decomposed in several spectral bands. As a consequence, after processing each buffer, the number of useful subbands is given by nb=2j-1, and processing m buffers we can form the pseudo-energy matrix E according  to (17), whose dimension is nb×m:
(17)E=[e1,1e1,2⋯e1,me2,1e2,2⋯e2,m⋮⋮⋯⋮enb,1enb,2⋯enb,m].From wavelet theory, the detailed coefficients are like noise with zero mean, and perturbations with rich spectral content will introduce the variations in the subbands pseudo-energy magnitude. So, the SVD/PCA analysis can be applied to obtain the hidden structure in the matrixE as established in (9).Note 1.
The idea behind this is to consider the noise and other noise-like variations as hidden structures of the sensors or, in other words, a sensor pattern.By decomposingE in to singular values (S), right eigenvectors (U), and left eigenvectors (V) and using (9),
(18)E=USVTand selecting its principal components (PCA) as stated in (10)–(13), (18) becomes
(19)E~=U~S~V~T,whereU~,S~, and V~T are partial matrices obtained from U,S, and VT according to PCA. The matrices E~ and E possess the same dimensions, but with different content.From PCA theory, it is known thatU is an orthonormal matrix that performs a rotation; then there exist, a null space projection matrix P, so that
(20)PU~=0.In the sequence, applyingP in (19), we obtain,
(21)PE~=PU~S~V~T,PE~=0.The meaning of (21) is that pseudo-energy vector obtained from buffer with certain number of samples at any time can be projected into a null space. So, the idea behind this is that the matrix P represents a signature of the sensor in the normal operational condition.Performing the operation given in (21) over m pseudo-energy vectors (see (17)) a dispersion about zero can be obtained (see 2D example in Figure 4). Now, assuming the rank of P as p, the Euclidean distance (dip) of the ith projection point is given by
(22)pdi=(xi1)2+(xi2)2+(xi3)2+⋯+(xip)2.Figure 4
2D example of the null space projection distances of the pseudo-energy vectors.As a consequence, the mean distance and its variance are given by(23)d-p=1m∑i=1mdip,(24)var(dp)=E[(dp-d-p)2],respectively.


#### 3.2.2. Noise Effect Reduction

The noise influence can be reduced by using a normalization technique as presented by de Oliveira [12]. In this technique the following model is assumed:
(25)y=cf+η,wherecf is a coefficient vector associated to the fault and η is the sensor noise vector.Assuming the average pseudo-energy at each subband as the autocorrelation ofy at τ=0 and that cf and η are uncorrelated, the energy equation becomes
(26)ey=ecf+eη,and the normalization is given by(27)e˘y=eyeη=ecfeη+1.Considering the variation of the noise energy between buffers, (27) can be rewritten as
(28)e˘y=eye-η≈ecfe-η+1,where thekth subband average energy (e-η) considering m buffers is given by
(29)(e-η)k=1m∑lm‍(eη)k,l,k=1:nbandl=1:m.Therefore, the normalized pseudo-energy matrix becomes(30)E=[e˘1,1e˘1,2⋯e˘1,me˘2,1e˘2,2⋯e˘2,m⋮⋮⋯⋮e˘nb,1e˘nb,2⋯e˘nb,m]and can be used instead ofE ((17)).


## 3.1. Fault Detection Algorithm

The detection of faults in an IMU with redundant sensors can be performed by PV analysis as defined in (3). 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 [13]. In the FD block, the recursive MF is used in order to remove impulsive noise and to improve the detection performed by χ2-CUSUM algorithm. In Figure 3 the FD processing structure in the left branch is shown. Under minimal redundancy configuration (four gyros), the FD algorithm can only indicate the fault occurrence. So, this information is used to start the identification process with high confidence of the fault occurrence. After the signal from gyros is processed by PV block according (3), the measurements are filtered by a second recursive MF (MF-2) to remove the spikes (or outliers) generated in the PV processing [6]. The filtered PV is applied to χ2-CUSUM block that performs the fault detection according to the algorithm and calibration presented in previous work [6].


## 3.2. Fault Isolation Algorithm

The fault isolation (FI) algorithm begins processing the data information stored in the buffers when the fault occurrence is indicated by FD block. The four buffers (one for each gyro) operate in FIFO mode. The operation of the FI algorithm is based on the analysis of the principal components in the energy subbands of the data stored in the buffers. The energy subbands are obtained by using the wavelet packet decomposition up to certain level. The idea behind this approach is to obtain asignature of the sensors in normal conditions given in terms of null space projection matrix and then use this matrix as kernel of FI analysis as demonstrated in the sequence.


### 3.2.1. FI Algorithm Model

The modeling of the FI algorithm begins with WP decomposition of a large amount of data from sensors divided in to pieces with properly defined size. The WP decomposition follows the schema illustrated in Figure2. Once the appropriated level of decomposition is reached, the respective subband elements are processed according to (16) that represents the average pseudo-energy level at each subband:
(16)ey=1N(ysb)T(ysb),whereN and ysb are the amounts of wavelet coefficients and respective wavelet coefficients vector at given subband.At each decomposition levelj, we have 2j subbands. The first subband (superscripts a,aa,aaa,…) is composed by approximation coefficients and represents the sensor state information (low-pass filter output), and this information will not be used in FI processing. The remaining subbands will form the pseudo-energy matrix (j≥2) because those subbands contain the noise information decomposed in several spectral bands. As a consequence, after processing each buffer, the number of useful subbands is given by nb=2j-1, and processing m buffers we can form the pseudo-energy matrix E according  to (17), whose dimension is nb×m:
(17)E=[e1,1e1,2⋯e1,me2,1e2,2⋯e2,m⋮⋮⋯⋮enb,1enb,2⋯enb,m].From wavelet theory, the detailed coefficients are like noise with zero mean, and perturbations with rich spectral content will introduce the variations in the subbands pseudo-energy magnitude. So, the SVD/PCA analysis can be applied to obtain the hidden structure in the matrixE as established in (9).Note 1.
The idea behind this is to consider the noise and other noise-like variations as hidden structures of the sensors or, in other words, a sensor pattern.By decomposingE in to singular values (S), right eigenvectors (U), and left eigenvectors (V) and using (9),
(18)E=USVTand selecting its principal components (PCA) as stated in (10)–(13), (18) becomes
(19)E~=U~S~V~T,whereU~,S~, and V~T are partial matrices obtained from U,S, and VT according to PCA. The matrices E~ and E possess the same dimensions, but with different content.From PCA theory, it is known thatU is an orthonormal matrix that performs a rotation; then there exist, a null space projection matrix P, so that
(20)PU~=0.In the sequence, applyingP in (19), we obtain,
(21)PE~=PU~S~V~T,PE~=0.The meaning of (21) is that pseudo-energy vector obtained from buffer with certain number of samples at any time can be projected into a null space. So, the idea behind this is that the matrix P represents a signature of the sensor in the normal operational condition.Performing the operation given in (21) over m pseudo-energy vectors (see (17)) a dispersion about zero can be obtained (see 2D example in Figure 4). Now, assuming the rank of P as p, the Euclidean distance (dip) of the ith projection point is given by
(22)pdi=(xi1)2+(xi2)2+(xi3)2+⋯+(xip)2.Figure 4
2D example of the null space projection distances of the pseudo-energy vectors.As a consequence, the mean distance and its variance are given by(23)d-p=1m∑i=1mdip,(24)var(dp)=E[(dp-d-p)2],respectively.


### 3.2.2. Noise Effect Reduction

The noise influence can be reduced by using a normalization technique as presented by de Oliveira [12]. In this technique the following model is assumed:
(25)y=cf+η,wherecf is a coefficient vector associated to the fault and η is the sensor noise vector.Assuming the average pseudo-energy at each subband as the autocorrelation ofy at τ=0 and that cf and η are uncorrelated, the energy equation becomes
(26)ey=ecf+eη,and the normalization is given by(27)e˘y=eyeη=ecfeη+1.Considering the variation of the noise energy between buffers, (27) can be rewritten as
(28)e˘y=eye-η≈ecfe-η+1,where thekth subband average energy (e-η) considering m buffers is given by
(29)(e-η)k=1m∑lm‍(eη)k,l,k=1:nbandl=1:m.Therefore, the normalized pseudo-energy matrix becomes(30)E=[e˘1,1e˘1,2⋯e˘1,me˘2,1e˘2,2⋯e˘2,m⋮⋮⋯⋮e˘nb,1e˘nb,2⋯e˘nb,m]and can be used instead ofE ((17)).


## 3.2.1. FI Algorithm Model

The modeling of the FI algorithm begins with WP decomposition of a large amount of data from sensors divided in to pieces with properly defined size. The WP decomposition follows the schema illustrated in Figure2. Once the appropriated level of decomposition is reached, the respective subband elements are processed according to (16) that represents the average pseudo-energy level at each subband:
(16)ey=1N(ysb)T(ysb),whereN and ysb are the amounts of wavelet coefficients and respective wavelet coefficients vector at given subband.At each decomposition levelj, we have 2j subbands. The first subband (superscripts a,aa,aaa,…) is composed by approximation coefficients and represents the sensor state information (low-pass filter output), and this information will not be used in FI processing. The remaining subbands will form the pseudo-energy matrix (j≥2) because those subbands contain the noise information decomposed in several spectral bands. As a consequence, after processing each buffer, the number of useful subbands is given by nb=2j-1, and processing m buffers we can form the pseudo-energy matrix E according  to (17), whose dimension is nb×m:
(17)E=[e1,1e1,2⋯e1,me2,1e2,2⋯e2,m⋮⋮⋯⋮enb,1enb,2⋯enb,m].From wavelet theory, the detailed coefficients are like noise with zero mean, and perturbations with rich spectral content will introduce the variations in the subbands pseudo-energy magnitude. So, the SVD/PCA analysis can be applied to obtain the hidden structure in the matrixE as established in (9).Note 1.
The idea behind this is to consider the noise and other noise-like variations as hidden structures of the sensors or, in other words, a sensor pattern.By decomposingE in to singular values (S), right eigenvectors (U), and left eigenvectors (V) and using (9),
(18)E=USVTand selecting its principal components (PCA) as stated in (10)–(13), (18) becomes
(19)E~=U~S~V~T,whereU~,S~, and V~T are partial matrices obtained from U,S, and VT according to PCA. The matrices E~ and E possess the same dimensions, but with different content.From PCA theory, it is known thatU is an orthonormal matrix that performs a rotation; then there exist, a null space projection matrix P, so that
(20)PU~=0.In the sequence, applyingP in (19), we obtain,
(21)PE~=PU~S~V~T,PE~=0.The meaning of (21) is that pseudo-energy vector obtained from buffer with certain number of samples at any time can be projected into a null space. So, the idea behind this is that the matrix P represents a signature of the sensor in the normal operational condition.Performing the operation given in (21) over m pseudo-energy vectors (see (17)) a dispersion about zero can be obtained (see 2D example in Figure 4). Now, assuming the rank of P as p, the Euclidean distance (dip) of the ith projection point is given by
(22)pdi=(xi1)2+(xi2)2+(xi3)2+⋯+(xip)2.Figure 4
2D example of the null space projection distances of the pseudo-energy vectors.As a consequence, the mean distance and its variance are given by(23)d-p=1m∑i=1mdip,(24)var(dp)=E[(dp-d-p)2],respectively.


## 3.2.2. Noise Effect Reduction

The noise influence can be reduced by using a normalization technique as presented by de Oliveira [12]. In this technique the following model is assumed:
(25)y=cf+η,wherecf is a coefficient vector associated to the fault and η is the sensor noise vector.Assuming the average pseudo-energy at each subband as the autocorrelation ofy at τ=0 and that cf and η are uncorrelated, the energy equation becomes
(26)ey=ecf+eη,and the normalization is given by(27)e˘y=eyeη=ecfeη+1.Considering the variation of the noise energy between buffers, (27) can be rewritten as
(28)e˘y=eye-η≈ecfe-η+1,where thekth subband average energy (e-η) considering m buffers is given by
(29)(e-η)k=1m∑lm‍(eη)k,l,k=1:nbandl=1:m.Therefore, the normalized pseudo-energy matrix becomes(30)E=[e˘1,1e˘1,2⋯e˘1,me˘2,1e˘2,2⋯e˘2,m⋮⋮⋯⋮e˘nb,1e˘nb,2⋯e˘nb,m]and can be used instead ofE ((17)).


## 4. Results

### 4.1. FI Algorithm

The analyses of algorithm performance were made upon two different data sets: One for analysis and calibration (alpha) and the other for validation (beta). Using the model developed in Section3.2, wavelet db4 (Daubechies 4), three decomposition levels (j=3), and a time series (alpha) with 51200 samples divided in to 200 buffers of dimension 256 samples, the matrices E7×200 and E˘7×200 were constructed. The dimension of the matrices is given by number of subbands (nb=2j-1) and the number of buffers. In this case, the rank of the matrices is seven (seven singular values). From those matrices with seven singular values four components were identified that represent more than 90% of the total energy (four principal components). With rank seven and using the first four greatest singular values, the matrix P (or P˘) will present the dimension 3×7 (see (A.5)–(A.8)), and consequently, a tridimensional projection (x1,x,2 and x3) according to the idea represented in Figure 4. In order to visualize, efficiency of the proposed model were performed the following operations: (a)
rotation:(U~)TE~,(b)
normalized rotation:(U˘~)TE˘~,(c)
parity space projection:PE~,(d)
normalized parity space projection:P˘E˘~,whose tridimensional projections are shown in Figure5 as scatter plot. In those four projections, the increase of the algorithm performance is straightforward  when the parity space projections approach is used (Figures 4(c) and 4(d)) compared to rotations (Figures 4(a) and 4(b)). Using (24), a 4σ spherical threshold was established in order to verify the consistency of the proposed model in terms of probability of false alarm (PFA), probability of miss detection (PMD), and correct decision (CD). In Figure 6 the tridimensional projection of the normalized parity space of the sensor noise decomposed by wavelet packet db4 at level j=3. As expected, the most of noise projections (about 97%) remains inside of the sphere, indicating a low occurence of false alarm. With a step bias fault randomly inserted in time into all of 200 buffers, it is expected the indication of 200 fault projections. In Figure 7 it is shown this operation comparing the noise and fault projections for the case of normalized parity space. The fault is of step bias kind with magnitude 0.8°/s.Scatter plot of 200 projections. (a) Rotation(U~)TE~; (b) normalized rotation (U˘~)TE˘~; (c) Parity space projection PE~; and (d) normalized Parity space projection P˘E˘~.
(a)
(b)
(c)
(d)Figure 6
Normalized parity space tridimensional projection of the sensor (gyro #2) noise decomposed by wavelet packet db4 at levelj=3 with 4σ spherical threshold.Figure 7
Normalized parity space tridimensional projection of the sensor (gyro#2) noise (black) and step bias fault (red) decomposed by wavelet packet db4 at level j=3 with 4σ spherical threshold.After extensive analyses comparing wavelet families (db2, db4, db6, haar, symlet 3, coiflet 1, and biorthogonal 3.1), decomposition level, buffer dimension, and fault level, it is summarized in Figures8 and 9 the best set of parameters to be used in the construction of the parity space projection matrix. Therefore, the best parameters found in order to used in the validation tests are WP db6, buffer size 64, and decomposition level j=3. The respective projection matrices (P˘) to be applied in validation tests are given by (A.5)–(A.8). So, using these parameters in another data set obtained with different movements of the IMU, it was found the results summarized in Table 2, where the low performance of gyro #4 is explained as a result of its higher noise level with respect  to the others sensors.Table 2
Validation results using the parameters obtained in the calibration over a different data set. Parameters: buffer = 64; WP db6;j=3; step bias fault = 0,6°/s; Pˇ given by (A.5)–(A.8).
gyro
%
Run
Mean
r1
r2
r3
r4
r5
r6
r7
r8
r9
r10
#1
PFA
2,6316
2,6316
2,6316
2,6316
2,6316
2,6316
2,6316
2,6316
2,6316
2,6316
2,6316
PMD
0,2506
0,2506
0,2506
0,2506
0,2506
0,5013
0,1253
0
0,5013
0,3759
0,2757
CD
99,7494
99,7494
99,7494
99,7494
99,7494
99,4987
99,8747
100
99,4987
99,6241
99,7243 ± 0,1541
#2
PFA
1,8797
1,8797
1,8797
1,8797
1,8797
1,8797
1,8797
1,8797
1,8797
1,8797
1,8797
PMD
0,3760
0,5013
0,5013
0,6266
0,6266
0,3760
0,5013
0,6266
0,6266
0,8772
0,5640
CD
99,6241
99,4987
99,4987
99,3734
99,3734
99,6241
99,4987
99,3734
99,3734
99,1228
99,4361 ± 0,1477
#3
PFA
2,8822
2,8822
2,8822
2,8822
2,8822
2,8822
2,8822
2,8822
2,8822
2,8822
2,8822
PMD
0,2506
0
0
0
0,5013
0,2506
0
0
0,6266
0,1253
0,1754
CD
99,7494
100
100
100
99,4987
99,7494
100
100
99,3734
99,8747
99,8246 ± 0,2303
#4
PFA
4,01
4,01
4,01
4,01
4,01
4,01
4,01
4,01
4,01
4,01
4,01
PMD
6,3910
6,1404
4,8872
5,3883
5,5138
4,7619
5,8897
4,8872
6,1404
4,3860
5,43859
CD
93,6090
93,8596
95,1128
94,6115
94,4862
95,2381
94,1103
95,1128
93,8596
95,6140
94,5614 ± 0,6894Figure 8
PFA, PMD and CD results for seven different wavelet families applied to gyro#2 data, considering the following parameters: buffer=256; j=3; step fault=0.8°/s.Figure 9
PFA, PMD and CD results for buffers dimensions 64, 128 and 256 applied to gyro#4 data, considering the following parameters: waveler db6; j=3; step fault=0.8°/s.


## 4.1. FI Algorithm

The analyses of algorithm performance were made upon two different data sets: One for analysis and calibration (alpha) and the other for validation (beta). Using the model developed in Section3.2, wavelet db4 (Daubechies 4), three decomposition levels (j=3), and a time series (alpha) with 51200 samples divided in to 200 buffers of dimension 256 samples, the matrices E7×200 and E˘7×200 were constructed. The dimension of the matrices is given by number of subbands (nb=2j-1) and the number of buffers. In this case, the rank of the matrices is seven (seven singular values). From those matrices with seven singular values four components were identified that represent more than 90% of the total energy (four principal components). With rank seven and using the first four greatest singular values, the matrix P (or P˘) will present the dimension 3×7 (see (A.5)–(A.8)), and consequently, a tridimensional projection (x1,x,2 and x3) according to the idea represented in Figure 4. In order to visualize, efficiency of the proposed model were performed the following operations: (a)
rotation:(U~)TE~,(b)
normalized rotation:(U˘~)TE˘~,(c)
parity space projection:PE~,(d)
normalized parity space projection:P˘E˘~,whose tridimensional projections are shown in Figure5 as scatter plot. In those four projections, the increase of the algorithm performance is straightforward  when the parity space projections approach is used (Figures 4(c) and 4(d)) compared to rotations (Figures 4(a) and 4(b)). Using (24), a 4σ spherical threshold was established in order to verify the consistency of the proposed model in terms of probability of false alarm (PFA), probability of miss detection (PMD), and correct decision (CD). In Figure 6 the tridimensional projection of the normalized parity space of the sensor noise decomposed by wavelet packet db4 at level j=3. As expected, the most of noise projections (about 97%) remains inside of the sphere, indicating a low occurence of false alarm. With a step bias fault randomly inserted in time into all of 200 buffers, it is expected the indication of 200 fault projections. In Figure 7 it is shown this operation comparing the noise and fault projections for the case of normalized parity space. The fault is of step bias kind with magnitude 0.8°/s.Scatter plot of 200 projections. (a) Rotation(U~)TE~; (b) normalized rotation (U˘~)TE˘~; (c) Parity space projection PE~; and (d) normalized Parity space projection P˘E˘~.
(a)
(b)
(c)
(d)Figure 6
Normalized parity space tridimensional projection of the sensor (gyro #2) noise decomposed by wavelet packet db4 at levelj=3 with 4σ spherical threshold.Figure 7
Normalized parity space tridimensional projection of the sensor (gyro#2) noise (black) and step bias fault (red) decomposed by wavelet packet db4 at level j=3 with 4σ spherical threshold.After extensive analyses comparing wavelet families (db2, db4, db6, haar, symlet 3, coiflet 1, and biorthogonal 3.1), decomposition level, buffer dimension, and fault level, it is summarized in Figures8 and 9 the best set of parameters to be used in the construction of the parity space projection matrix. Therefore, the best parameters found in order to used in the validation tests are WP db6, buffer size 64, and decomposition level j=3. The respective projection matrices (P˘) to be applied in validation tests are given by (A.5)–(A.8). So, using these parameters in another data set obtained with different movements of the IMU, it was found the results summarized in Table 2, where the low performance of gyro #4 is explained as a result of its higher noise level with respect  to the others sensors.Table 2
Validation results using the parameters obtained in the calibration over a different data set. Parameters: buffer = 64; WP db6;j=3; step bias fault = 0,6°/s; Pˇ given by (A.5)–(A.8).
gyro
%
Run
Mean
r1
r2
r3
r4
r5
r6
r7
r8
r9
r10
#1
PFA
2,6316
2,6316
2,6316
2,6316
2,6316
2,6316
2,6316
2,6316
2,6316
2,6316
2,6316
PMD
0,2506
0,2506
0,2506
0,2506
0,2506
0,5013
0,1253
0
0,5013
0,3759
0,2757
CD
99,7494
99,7494
99,7494
99,7494
99,7494
99,4987
99,8747
100
99,4987
99,6241
99,7243 ± 0,1541
#2
PFA
1,8797
1,8797
1,8797
1,8797
1,8797
1,8797
1,8797
1,8797
1,8797
1,8797
1,8797
PMD
0,3760
0,5013
0,5013
0,6266
0,6266
0,3760
0,5013
0,6266
0,6266
0,8772
0,5640
CD
99,6241
99,4987
99,4987
99,3734
99,3734
99,6241
99,4987
99,3734
99,3734
99,1228
99,4361 ± 0,1477
#3
PFA
2,8822
2,8822
2,8822
2,8822
2,8822
2,8822
2,8822
2,8822
2,8822
2,8822
2,8822
PMD
0,2506
0
0
0
0,5013
0,2506
0
0
0,6266
0,1253
0,1754
CD
99,7494
100
100
100
99,4987
99,7494
100
100
99,3734
99,8747
99,8246 ± 0,2303
#4
PFA
4,01
4,01
4,01
4,01
4,01
4,01
4,01
4,01
4,01
4,01
4,01
PMD
6,3910
6,1404
4,8872
5,3883
5,5138
4,7619
5,8897
4,8872
6,1404
4,3860
5,43859
CD
93,6090
93,8596
95,1128
94,6115
94,4862
95,2381
94,1103
95,1128
93,8596
95,6140
94,5614 ± 0,6894Figure 8
PFA, PMD and CD results for seven different wavelet families applied to gyro#2 data, considering the following parameters: buffer=256; j=3; step fault=0.8°/s.Figure 9
PFA, PMD and CD results for buffers dimensions 64, 128 and 256 applied to gyro#4 data, considering the following parameters: waveler db6; j=3; step fault=0.8°/s.


## 5. Conclusion

In this paper, it was developed a FDI algorithm based on signal processing techniques applied to IMU with minimal redundancy of fiber optic gyros. This unity was built with low-quality FOGs and a prototype of acquisition system that generates a high quantity of impulsive noise, justifying the use of recursive median filter. The first of median filters (MF-1) was applied to remove the high-level impulses, outliers with huge magnitude arising from data acquisition process. This filter is adjusted to introduce a short delay. The second filter (MF-2) was used to reduce the low-level impulses arising from parity vector operations and,consequently, to improve the fault detection performance. However, this filter adds a new delay in the detection process, being greater than MF-1. The fault detection part of the algorithm is used to guarantee the reliability of the isolation part, avoiding the actions when the fault occurrence is not confirmed byχ2-CUSUM algorithm. In the fault detection part a calibration referenced on fault isolation characteristics was developed, mainly the dimension of the buffer used in the wavelet decomposition. The fault detection results demonstrated the effectiveness of the algorithm at low level of fault with short delay (about 10 samples for 0.3°/s step bias fault) [6]. The fault isolation algorithm, based on parity space projection of subband energy obtained from wavelet packet decomposition, provided expressive results in terms of correct decisions (CD >94% for 0.6°/s step bias fault). This parity space projection of subband energy is an extension of clustering analysis based on PCA. The step bias fault was used in both algorithms in order to verify the sensibility, but this algorithm can detect other kinds of faults that perform variations in the spectral contents of the sensors.



---

*Source: 869293-2013-11-21.xml*