# Research on Wavelet Singularity Detection Based Fault-Tolerant Federated Filtering Algorithm for INS/GPS/DVL Integrated Navigation System

Authors: Jianhua Cheng; Daidai Chen; Rene Landry Jr.; Lin Zhao; Dongxue Guan
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

Soft faults in navigation sensors will lead to the degradation of the accuracy and reliability of integrated navigation system. To solve this problem, a wavelet analysis and signal singularities based soft fault detection method are given out. To find signal singularities and detect the faults, the modulus maxima values are calculated after the wavelet transform of original signal. By calculating the Lipschitz exponent using the modulus maxima value at the fault point, the fault types are distinguished. Then, a fault-tolerant federated filtering algorithm for the calibration of INS/GPS/DVL integrated navigation system is proposed. Simulations are conducted and results show that sensor soft faults can be detected accurately. By effectively isolating the fault and refactoring information, the accuracy and reliability of navigation system are improved.


---

## Body

## 1. Introduction

Higher standards for the accuracy and reliability of navigation systems are required with the development of modern ship, aircraft, and land vehicle. Device redundancy and estimator based analytical redundancy are the two most effective methods of improving the accuracy and reliability of system [1].Device redundancy appears to be a simple and efficient approach to improve the performance of navigation system [2]. Both system level and sensor level redundancy are considerable, and most of the system redundancies are in the form of integrated navigation system. And INS/GPS/DVL integrated navigation system is widely used in vehicles. Inertial navigation system (INS) has been usually selected as the core of integrated systems for its outstanding characters, such as high-accuracy, self-contained, and jam-proof. Meanwhile, global position system (GPS) and Doppler velocity log (DVL) commonly play the external auxiliary roles, because their navigation errors do not accumulate with time [3]. They also provide the position and velocity for the calibration and reset of INS. That means that the performance of GPS and DVL will affect the calibration accuracy of INS directly as well as the entire INS/GPS/DVL integrated navigation system. As for the entire calibration process for INS, it is necessary to guarantee the accuracy and reliability of external information [4]. Otherwise, the gyroscope drifts cannot be estimated and the error of external information will be introduced.Estimator based analytical redundancy method is another effective solution [5, 6]. The estimated information predicted by using system mathematical model is introduced, then conduct data fusion by estimator (i.e., Kalman filter) so that almost all outputs of different kinds of sensors can be fully used. However, the fault of sensors might happen, which must be detected and processed. Kalman filter based fault detection methods are commonly used in integrated navigation system, especially χ
2 detection method.Inχ
2 detection method, faults are detected by comparing the assumptive mean and variance with those of constructedn-dimension Gaussian variables, which is easy to realize. However, this method insensitive and large-delay for the detection of soft faults, and the type of fault cannot be accurately determined. In the integrated system, there are two most common sensor faults: (
1
) step error caused by clock jumping and power fluctuation in GPS; (
2
) slope error caused by clock drift and the inaccuracy of rail parameter model [7].To solve the problems above, a sensor fault detection method based on wavelet analysis and signal singularities is put forward, and a fault-tolerant federated filter for INS/GPS/DVL integrated navigation system is designed. In Section2, signal fault detection method based on wavelet modulus maxima is introduced, and the relationships between signal singularity and modulus maxima are given out, and by searching modulus maxima the fault can be accurately positioned. In Section 3, fault type detection method based on the calculation of signal Lipschitz exponent is proposed. In Section 4, by utilizing the proposed fault detection method, a fault-tolerant federated filter for the calibration progress in INS/GPS/DVL integrated navigation system is designed. And in Section 5, simulations for fault-tolerant federated filter in calibration application are conducted. Finally, the conclusions are given in Section 6.


## 2. Signal Singularity Detection Method Based on Wavelet Analysis

### 2.1. Wavelet Transform and Modulus Maxima

Wavelet analysis is an efficient signal analysis tool and has been widely used in signal process [8, 9]. Wavelet transform appears to be extremely sensitive for the signal singularities, which makes it highly appropriate for the fault signal detection. In this transform, the original signal is processed by using wavelet as follows:
(1)
ψ
s
,
u
(
t
)
=
1
|
s
|
ψ
(
t
-
u
s
)
s
,
u
∈
R
,
s
≠
0
,
where the ψ
(
t
) is mother wavelet, s is scale parameter, and u is translation parameter.By (1), with the decrease of |
s
|, the locality and resolution in time domain of wavelet function are improved, and resolution in frequency domain is reduced.Assuming that waveletψ
(
t
) and the sensor output signal f
(
t
)
∈
L
2
(
R
) are both real functions, scale parameter s
>
0. The continuous wavelet transform for f
(
t
) can be calculated as follows:
(2)
W
f
(
s
,
u
)
=
s
1
/
2
(
f
*
θ
-
s
)
(
u
)
=
s
1
/
2
d
d
u
(
f
*
θ
-
s
)
(
u
)
,
where “*” denotes the convolution operation and θ
-
s is the complex conjugate of θ
s.As known from (3), under the determined scale s
0, we could obtain the local maxima of W
f
(
s
,
u
) when the derivative ∂
W
f
(
s
0
,
u
)
/
∂
u
=
0 at u
=
u
0. If any points u within the neighbourhood of u
0 satisfy equation W
f
(
s
0
,
u
)
≤
W
f
(
s
0
,
u
0
) and any points u within the left or right neighbourhood of u
0 satisfy equation W
f
(
s
0
,
u
)
<
W
f
(
s
0
,
u
0
), we call the point (
s
0
,
u
0
) as the extremum of the wavelet transform modulus W
f
(
s
,
u
) on the scale s
0 and |
W
f
(
s
0
,
u
0
)
| as the modulus maxima for the wavelet transform for the original f
(
t
).Due to the high-sampling ratio and nature discrete characters for the output signals of navigation systems, the continuous wavelet transform leads to extremely heavy calculation burden. It seems impossible to implement continuous wavelet transform in practical signal processing. Therefore, binary wavelet transform is selected in this paper, which can be seen as a compromise between continuous and discrete wavelet transforms. And lets
=
2
j; the binary wavelet can be expressed as
(3)
ψ
2
j
,
u
(
t
)
=
2
-
j
/
2
ψ
[
2
-
j
(
t
-
u
)
]
j
∈
Z
,
u
∈
R
.Accordingly, the binary wavelet transform is(4)
W
f
(
2
j
,
u
)
=
2
-
j
/
2
∫
-
∞
+
∞
f
(
t
)
ψ
[
2
-
j
(
t
-
u
s
)
]
d
t
.Signal singularities, which are the most important local information of signals, can be inferred by searching the modulus maxima. By calculating the modulus maxima on different scales, the entire modulus maxima curve can be obtained, and then we can observe the local signal singularities. The relationships among the original signal and the modulus maxima are shown in Figure1.Figure 1
Signal singularities detection based on wavelet modulus maxima.As shown in Figure1, if the characteristics of signal changes, their corresponding wavelet modulus appears as maxima. This character can be utilized to detect and position the sensor faults.


### 2.2. Fault Positioning Based on Modulus Maxima Searching

To detect the singularity of fault signal, we have to obtain the curves of modulus maxima. And ad hoc algorithm is commonly used in searching the points of modulus maxima [10].The specific steps of the algorithm are shown as follows.(1)
Choose Morlet wavelet basis and scale2
j and process the sensor output signal by using binary discrete wavelet transform. The chosen scale should guarantee that the signal extrema are dominated on the largest scale and keep all the important singularities of signal. In order to make the output signal meet the requirement of real time for the integrated navigation system, j
=
5 is set.
(2)
Assuming that the largest amplitude of extrema isM on the largest scale, the extrema whose amplitudes are smaller than M
/
j should be eliminated, that is, because these extrema are mainly caused by noise.
(3)
For arbitrary extremumx
0 on scale 2
j, its neighbour points are x
1 and x
2, and the manifold-points of x
1 are x
1
′ should be in the interval L
=
[
max
⁡
(
x
1
,
x
1
′
)
,
x
2
]. Among the points (
a
1
,
a
2
,
…
,
a
n
) have the same sign with x
0 in the interval L, if a
k satisfies |
a
k
-
x
0
|
≤
|
a
i
-
x
0
|
/
3, i
=
1,2
,
…
,
n
,
k
≠
i, a
k is the manifold-point of x
0. Otherwise, the point with the largest amplitude is the manifold point of x
0. However, once the amplitude of the manifold point is twice that of x
0, this point will be treated as the extrema of noise and have to be eliminated [11].
(4)
Repeat the process above until the scale is2
2.
(5)
Eliminate the maximum points on scale2
1, and use the extrema on scale 2
2. So far, the search for the modulus maxima on scales 2
j
~
2
1 is completed.Through the steps above, the modulus maxima on each scale can be obtained as well as the modulus maxima curve. According to the modulus maxima on fine scale, we can accurately position the singular points and the moments of sensor faults that occurred. However, only with these wavelet modulus maxima, we still could not identify the types of sensor faults.


## 2.1. Wavelet Transform and Modulus Maxima

Wavelet analysis is an efficient signal analysis tool and has been widely used in signal process [8, 9]. Wavelet transform appears to be extremely sensitive for the signal singularities, which makes it highly appropriate for the fault signal detection. In this transform, the original signal is processed by using wavelet as follows:
(1)
ψ
s
,
u
(
t
)
=
1
|
s
|
ψ
(
t
-
u
s
)
s
,
u
∈
R
,
s
≠
0
,
where the ψ
(
t
) is mother wavelet, s is scale parameter, and u is translation parameter.By (1), with the decrease of |
s
|, the locality and resolution in time domain of wavelet function are improved, and resolution in frequency domain is reduced.Assuming that waveletψ
(
t
) and the sensor output signal f
(
t
)
∈
L
2
(
R
) are both real functions, scale parameter s
>
0. The continuous wavelet transform for f
(
t
) can be calculated as follows:
(2)
W
f
(
s
,
u
)
=
s
1
/
2
(
f
*
θ
-
s
)
(
u
)
=
s
1
/
2
d
d
u
(
f
*
θ
-
s
)
(
u
)
,
where “*” denotes the convolution operation and θ
-
s is the complex conjugate of θ
s.As known from (3), under the determined scale s
0, we could obtain the local maxima of W
f
(
s
,
u
) when the derivative ∂
W
f
(
s
0
,
u
)
/
∂
u
=
0 at u
=
u
0. If any points u within the neighbourhood of u
0 satisfy equation W
f
(
s
0
,
u
)
≤
W
f
(
s
0
,
u
0
) and any points u within the left or right neighbourhood of u
0 satisfy equation W
f
(
s
0
,
u
)
<
W
f
(
s
0
,
u
0
), we call the point (
s
0
,
u
0
) as the extremum of the wavelet transform modulus W
f
(
s
,
u
) on the scale s
0 and |
W
f
(
s
0
,
u
0
)
| as the modulus maxima for the wavelet transform for the original f
(
t
).Due to the high-sampling ratio and nature discrete characters for the output signals of navigation systems, the continuous wavelet transform leads to extremely heavy calculation burden. It seems impossible to implement continuous wavelet transform in practical signal processing. Therefore, binary wavelet transform is selected in this paper, which can be seen as a compromise between continuous and discrete wavelet transforms. And lets
=
2
j; the binary wavelet can be expressed as
(3)
ψ
2
j
,
u
(
t
)
=
2
-
j
/
2
ψ
[
2
-
j
(
t
-
u
)
]
j
∈
Z
,
u
∈
R
.Accordingly, the binary wavelet transform is(4)
W
f
(
2
j
,
u
)
=
2
-
j
/
2
∫
-
∞
+
∞
f
(
t
)
ψ
[
2
-
j
(
t
-
u
s
)
]
d
t
.Signal singularities, which are the most important local information of signals, can be inferred by searching the modulus maxima. By calculating the modulus maxima on different scales, the entire modulus maxima curve can be obtained, and then we can observe the local signal singularities. The relationships among the original signal and the modulus maxima are shown in Figure1.Figure 1
Signal singularities detection based on wavelet modulus maxima.As shown in Figure1, if the characteristics of signal changes, their corresponding wavelet modulus appears as maxima. This character can be utilized to detect and position the sensor faults.


## 2.2. Fault Positioning Based on Modulus Maxima Searching

To detect the singularity of fault signal, we have to obtain the curves of modulus maxima. And ad hoc algorithm is commonly used in searching the points of modulus maxima [10].The specific steps of the algorithm are shown as follows.(1)
Choose Morlet wavelet basis and scale2
j and process the sensor output signal by using binary discrete wavelet transform. The chosen scale should guarantee that the signal extrema are dominated on the largest scale and keep all the important singularities of signal. In order to make the output signal meet the requirement of real time for the integrated navigation system, j
=
5 is set.
(2)
Assuming that the largest amplitude of extrema isM on the largest scale, the extrema whose amplitudes are smaller than M
/
j should be eliminated, that is, because these extrema are mainly caused by noise.
(3)
For arbitrary extremumx
0 on scale 2
j, its neighbour points are x
1 and x
2, and the manifold-points of x
1 are x
1
′ should be in the interval L
=
[
max
⁡
(
x
1
,
x
1
′
)
,
x
2
]. Among the points (
a
1
,
a
2
,
…
,
a
n
) have the same sign with x
0 in the interval L, if a
k satisfies |
a
k
-
x
0
|
≤
|
a
i
-
x
0
|
/
3, i
=
1,2
,
…
,
n
,
k
≠
i, a
k is the manifold-point of x
0. Otherwise, the point with the largest amplitude is the manifold point of x
0. However, once the amplitude of the manifold point is twice that of x
0, this point will be treated as the extrema of noise and have to be eliminated [11].
(4)
Repeat the process above until the scale is2
2.
(5)
Eliminate the maximum points on scale2
1, and use the extrema on scale 2
2. So far, the search for the modulus maxima on scales 2
j
~
2
1 is completed.Through the steps above, the modulus maxima on each scale can be obtained as well as the modulus maxima curve. According to the modulus maxima on fine scale, we can accurately position the singular points and the moments of sensor faults that occurred. However, only with these wavelet modulus maxima, we still could not identify the types of sensor faults.


## 3. Fault Type Detection Based on Lipschitz Exponent

### 3.1. Lipschitz Exponent Characteristics of Signal Faults

Singularity refers to the signal with discontinuous points or its derivative function is not continuous. Mathematically, signal singularity is described by Lipschitz exponent [12]. The type of sensor faults could be identified based on Lipschitz exponent.The Lipschitz exponentα of signal f
(
t
) at t
=
h
0 is defined as follows. There are two constants M, h
0
(
M
>
0
,
h
0
>
0
) and a polynomial g
n
(
h
). When h
<
h
0, inequation is established as
(5)
|
f
(
h
0
+
h
)
-
g
n
(
h
0
)
|
≤
M
|
h
|
α
.For any pointh
0, the polynomial g
n
(
h
) is uniquely determined. If f
(
t
) is n
(
n
=
⌊
α
⌋
) order continuously differentiable at point h
0 and its neighborhood, the g
n
(
h
) is equal to the Taylor expansion of f
(
t
) at point h
0. If function f
(
t
) has uniform Lipschitz exponent in any neighborhood of h
0 point and α
<
n, f
(
t
) is n order differentiable in this neighborhood of point h
0. Therefore, a great Lipschitz exponent α means better smoothness of signal and the smaller α characterizes the singularities of f
(
t
).As for the most typical sensor errors, the relationships between slope error, step error, random noise error, and Lipschitz exponents are given out.(1)
Assuming thatf
(
t
) is a step signal, we have
(6)
|
f
(
h
0
+
h
)
-
g
n
(
h
0
)
|
≤
O
|
h
|
0
.According to (5), the Lipschitz exponent α of step signal at h
0 point is 0.(2)
Assuming thatf
(
t
) is a slope signal, we have
(7)
|
f
(
h
0
+
h
)
-
g
n
(
h
0
)
|
≤
O
|
h
|
1
.According to (5), the Lipschitz exponent α of slope signal at h
0 point is 1.However, the Lipschitz exponents related to random noise are usually less than 0, because Gaussian white noise is a random distribution whose singularities are almost everywhere. On this condition, the Lipschitz exponents are negative and could be denoted asα
=
-
0.5
-
ε
,
∀
ε
>
0.


### 3.2. Calculation of Signal Lipschitz Exponent

By analyzing the relationship between the modulus maxima and signal singularity, it is found that the attenuation of|
W
f
(
s
,
u
)
| could be controlled by its local maxima. Therefore, according to the relationship between wavelet modulus maxima and signal Lipschitz exponent, the calculation of signal Lipschitz exponent can be realized [13].Assuming that the compactly supported set of wavelet functionψ is [
-
C
,
C
]
(
C
>
0
) and when s
<
s
0, all the modulus maxima points which converge to the point v are located in the conical surface, which is denoted as |
u
-
v
|
≤
C
s, as shown in Figure 2.Figure 2
The conical surface composed of modulus maxima points.Assuming the signalf
(
t
)
∈
L
2
(
R
) and its modulus maximum values distribute within a conical surface. That is
(8)
|
W
f
(
s
,
u
)
|
≤
A
a
α
+
0.5
·
|
1
+
|
t
-
t
0
s
|
α
|
,
where A is a constant. 								(9)
|
W
f
(
s
,
u
)
|
≤
A
a
α
+
0.5
.Calculating the logarithm for both sides of (9) and using binary wavelet transform, it can be obtained that
(10)
lo
g
2
|
W
f
(
2
j
,
t
j
)
|
≤
lo
g
2
A
+
j
(
α
+
0.5
)
.Equation (10) shows that the Lipschitz exponent value at point t
0 of f
(
t
) depends on the attenuation of |
W
f
(
s
,
u
)
| under fine scales. And when t
=
t
0, we get the following conclusions.(a)
If Lipschitz exponentα
>
0, the modulus maxima is proportional to the scale.
(b)
If Lipschitz exponentα
<
0, the modulus maxima is inversely proportional to the scale.
(c)
If Lipschitz exponentα
=
0, the modulus maxima is uncorrelated with the scale.And for the next scalej
+
1, we have
(11)
lo
g
2
|
W
f
(
2
j
+
1
,
t
j
+
1
)
|
≤
lo
g
2
A
+
(
j
+
1
)
(
α
+
0.5
)
.Assuming that the equalities hold up in both (10) and (11), the formula for the calculation of the Lipschitz exponents is as follows:
(12)
α
=
lo
g
2
|
W
f
(
2
j
+
1
,
t
j
+
1
)
|
-
lo
g
2
|
W
f
(
2
j
,
t
j
)
|
-
0.5
.Through the steps above, calculate the modulus maxima on every scale and draw the modulus maxima curve. According to the modulus maxima on fine scale, we can find the singular point position, namely, the moment of fault occurrence.Select the modulus maxima point on every scale calculated in step (6) and the logarithm of decomposition scale as ordinate and abscissa, respectively. The slope of curve is estimated by linear least squares, and Lipschitz exponent can be obtained by subtracting 0.5 from the curve slope. According to the Lipschitz exponent, we can distinguish the fault type, and the fault information is provided for fault detection and isolation module to dispose of the sensors with faults.


## 3.1. Lipschitz Exponent Characteristics of Signal Faults

Singularity refers to the signal with discontinuous points or its derivative function is not continuous. Mathematically, signal singularity is described by Lipschitz exponent [12]. The type of sensor faults could be identified based on Lipschitz exponent.The Lipschitz exponentα of signal f
(
t
) at t
=
h
0 is defined as follows. There are two constants M, h
0
(
M
>
0
,
h
0
>
0
) and a polynomial g
n
(
h
). When h
<
h
0, inequation is established as
(5)
|
f
(
h
0
+
h
)
-
g
n
(
h
0
)
|
≤
M
|
h
|
α
.For any pointh
0, the polynomial g
n
(
h
) is uniquely determined. If f
(
t
) is n
(
n
=
⌊
α
⌋
) order continuously differentiable at point h
0 and its neighborhood, the g
n
(
h
) is equal to the Taylor expansion of f
(
t
) at point h
0. If function f
(
t
) has uniform Lipschitz exponent in any neighborhood of h
0 point and α
<
n, f
(
t
) is n order differentiable in this neighborhood of point h
0. Therefore, a great Lipschitz exponent α means better smoothness of signal and the smaller α characterizes the singularities of f
(
t
).As for the most typical sensor errors, the relationships between slope error, step error, random noise error, and Lipschitz exponents are given out.(1)
Assuming thatf
(
t
) is a step signal, we have
(6)
|
f
(
h
0
+
h
)
-
g
n
(
h
0
)
|
≤
O
|
h
|
0
.According to (5), the Lipschitz exponent α of step signal at h
0 point is 0.(2)
Assuming thatf
(
t
) is a slope signal, we have
(7)
|
f
(
h
0
+
h
)
-
g
n
(
h
0
)
|
≤
O
|
h
|
1
.According to (5), the Lipschitz exponent α of slope signal at h
0 point is 1.However, the Lipschitz exponents related to random noise are usually less than 0, because Gaussian white noise is a random distribution whose singularities are almost everywhere. On this condition, the Lipschitz exponents are negative and could be denoted asα
=
-
0.5
-
ε
,
∀
ε
>
0.


## 3.2. Calculation of Signal Lipschitz Exponent

By analyzing the relationship between the modulus maxima and signal singularity, it is found that the attenuation of|
W
f
(
s
,
u
)
| could be controlled by its local maxima. Therefore, according to the relationship between wavelet modulus maxima and signal Lipschitz exponent, the calculation of signal Lipschitz exponent can be realized [13].Assuming that the compactly supported set of wavelet functionψ is [
-
C
,
C
]
(
C
>
0
) and when s
<
s
0, all the modulus maxima points which converge to the point v are located in the conical surface, which is denoted as |
u
-
v
|
≤
C
s, as shown in Figure 2.Figure 2
The conical surface composed of modulus maxima points.Assuming the signalf
(
t
)
∈
L
2
(
R
) and its modulus maximum values distribute within a conical surface. That is
(8)
|
W
f
(
s
,
u
)
|
≤
A
a
α
+
0.5
·
|
1
+
|
t
-
t
0
s
|
α
|
,
where A is a constant. 								(9)
|
W
f
(
s
,
u
)
|
≤
A
a
α
+
0.5
.Calculating the logarithm for both sides of (9) and using binary wavelet transform, it can be obtained that
(10)
lo
g
2
|
W
f
(
2
j
,
t
j
)
|
≤
lo
g
2
A
+
j
(
α
+
0.5
)
.Equation (10) shows that the Lipschitz exponent value at point t
0 of f
(
t
) depends on the attenuation of |
W
f
(
s
,
u
)
| under fine scales. And when t
=
t
0, we get the following conclusions.(a)
If Lipschitz exponentα
>
0, the modulus maxima is proportional to the scale.
(b)
If Lipschitz exponentα
<
0, the modulus maxima is inversely proportional to the scale.
(c)
If Lipschitz exponentα
=
0, the modulus maxima is uncorrelated with the scale.And for the next scalej
+
1, we have
(11)
lo
g
2
|
W
f
(
2
j
+
1
,
t
j
+
1
)
|
≤
lo
g
2
A
+
(
j
+
1
)
(
α
+
0.5
)
.Assuming that the equalities hold up in both (10) and (11), the formula for the calculation of the Lipschitz exponents is as follows:
(12)
α
=
lo
g
2
|
W
f
(
2
j
+
1
,
t
j
+
1
)
|
-
lo
g
2
|
W
f
(
2
j
,
t
j
)
|
-
0.5
.Through the steps above, calculate the modulus maxima on every scale and draw the modulus maxima curve. According to the modulus maxima on fine scale, we can find the singular point position, namely, the moment of fault occurrence.Select the modulus maxima point on every scale calculated in step (6) and the logarithm of decomposition scale as ordinate and abscissa, respectively. The slope of curve is estimated by linear least squares, and Lipschitz exponent can be obtained by subtracting 0.5 from the curve slope. According to the Lipschitz exponent, we can distinguish the fault type, and the fault information is provided for fault detection and isolation module to dispose of the sensors with faults.


## 4. Design for Fault-Tolerant Filtering Scheme

To maintain the high accuracy of INS/GPS/DVL integrated navigation system for long term operation, comprehensive calibration is one of the key technologies [14]. External information of INS such as position provided by GPS is introduced to estimate the gyroscope drifts and reset the navigation information. Usually, this progress might last for hours. Moreover, during the entire calibration phase, the accuracy of external information must be maintained at a stable level. However, the GPS signals appear to be easily influenced by dynamic environment or human disturbances.In order to guarantee the accuracy and reliability of INS/GPS/DVL integrated navigation system, a fault detection method based fault-tolerant filtering scheme for the calibration phase is presented. In this scheme, decentralized filter is adopted; the fault detection method mentioned above is used to process the outputs of the two subintegrated systems INS/GPS and INS/DVL. And the diagram of this fault-tolerant filtering algorithm is given in Figure3.Figure 3
The fault-tolerant federated filter based on fault detection method for INS/GPS/DVL integrated navigation system.As shown in Figure3, position measurementZ
IG is the position difference between GPS and INS outputs, and velocity measurementZ
ID is the velocity difference between DVL and INS outputs. Then,Z
IGandZ
IDare processed through the fault detection method based on binary wavelet transform and Lipschitz exponent. Once faults are detected, the faults can be quickly and accurately positioned and then the type of fault would be obtained. Consequently, the fault sensor would be isolated, and the corresponding output would be reconstructed to form measurement informationZ
IGandZ
ID which are sent to subfilters. That means that the integrated navigation system exhibits better fault-tolerant ability. With more precise external navigation information, the integrated system could normally output precise navigation information and even sensor soft fault occurs.


### 4.1. State and Measurement Equations for Subfilter

Considering that the outputs from the sensors are different, we use indirect method to estimate the navigation parameters. The state vector of INS is selected as(13)
X
I
=
[
δ
φ
I
δ
λ
I
δ
V
x
δ
V
y
α
β
γ
ε
x
ε
y
ε
z
ε
r
x
ε
r
y
ε
r
z
]
T
,
where δ
φ,  δ
λ are the position errors of INS; δ
V
x, δ
V
y are the velocity errors of INS; α, β and γ are attitude errors of INS; ε
x, ε
y and ε
z are the gyro constant drifts; ε
r
x, ε
r
y and ε
r
z are the gyro random drifts.The state vector of GPS is selected as(14)
X
G
=
[
δ
φ
G
δ
λ
G
]
T
,
where δ
φ
G, δ
λ
G are the position error of GPS.Position integrated mode is used for INS/GPS subfilter, so the state equation and measurement equation are given by(15)
[
X
˙
I
X
˙
G
]
=
[
F
I
0
0
F
G
]
[
X
I
X
G
]
+
[
W
I
W
G
]
,
Z
I
G
=
[
φ
I
-
φ
G
λ
I
-
λ
G
]
=
H
I
G
[
X
I
X
G
]
+
V
I
G
,
where F
I is INS state transition matrix in [15]; W
I, W
G are the system noise matrixes for INS and GPS, respectively; V
I
G is measurement noise matrix; and H
I
G is state transition matrix for INS/GPS measurement equation. Consider
(16)
W
I
=
[
ω
g
x
ω
g
y
ω
g
z
]
T
,
W
G
=
[
ω
φ
ω
λ
]
T
,
where ω
g
x, ω
g
y and ω
g
z are system white noises of INS and ω
φ, ω
λ are system white noise of GPS. Consider
(17)
H
I
G
=
[
1
0
0
0
0
0
0
0
0
0
0
0
0
-
1
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
-
1
]
.Using first order Markov process to complete the approximate fitting, GPS state transition matrixF
G can be expressed as follows:
(18)
F
G
=
diag
⁡
(
1
τ
GPS
φ
,
1
τ
GPS
λ
)
,
where τ
GPS
φ and τ
GPS
λ are the correlation times.Velocity matching method is used in INS/DVL subfilter, and the state equation and measurement equation are given by(19)
[
X
˙
I
X
˙
D
]
=
[
F
I
0
0
F
D
]
[
X
I
X
D
]
+
[
W
I
W
D
]
,
Z
I
D
=
[
V
I
x
-
V
D
x
V
I
y
-
V
D
y
]
=
H
I
D
[
X
I
X
D
]
+
V
I
D
,
where X
D
=
[
δ
V
D
x
δ
V
D
y
δ
k
]
T is state vector of DVL velocity estimation error, δ
V
D
x, δ
V
D
y are velocity error, and δ
k is scale error. F
D is DVL state transition matrix:
(20)
F
D
=
diag
⁡
(
1
τ
D
x
,
1
τ
D
y
,
1
τ
k
)
,
where τ
D
x, τ
D
y and τ
k are Markov correlation time, W
D is system noise matrix, ω
D
x, ω
D
y are driven white noise, and V
I
D is measure noise matrix. Consider
(21)
W
G
=
[
ω
D
x
ω
D
y
]
T
,
H
I
D
=
[
0
0
1
0
0
0
-
V
y
0
0
0
0
0
0
-
1
0
-
V
x
0
0
0
1
0
0
V
x
0
0
0
0
0
0
0
-
1
-
V
y
]
.


### 4.2. INS/GPS/DVL Global Filter

To make the integrated filter realizable for computer calculation, it is necessary to discretize the state equation and measurement equation mentioned above and then we have(22)
X
k
=
Φ
k
,
k
-
1
X
k
+
Γ
k
-
1
W
k
-
1
,
Z
k
=
H
k
X
k
+
V
k
.Assuming that the local estimates are unrelated, the global estimate can be expressed as(23)
X
^
g
=
P
g
(
P
11
-
1
X
^
1
+
P
22
-
1
X
^
2
)
,
P
g
=
(
P
11
-
1
+
P
22
-
1
)
-
1
.Local estimates (X
^
1
,
X
^
2
) from the subfilters and their covariance matrix (P
1
,
P
2
) are sent to master filter, and the obtained results would integrate with master filter estimates to get the global optimal estimates. Besides, global estimate X
^
g  and its covariance matrix P
g, which are amplified to be β
i
-
1
P
g
(
β
i
≤
1
), are feed back to sub filters, and estimates of subfilter are reset as
(24)
X
^
i
=
X
^
g
,
P
i
i
=
β
i
-
1
P
g
.


## 4.1. State and Measurement Equations for Subfilter

Considering that the outputs from the sensors are different, we use indirect method to estimate the navigation parameters. The state vector of INS is selected as(13)
X
I
=
[
δ
φ
I
δ
λ
I
δ
V
x
δ
V
y
α
β
γ
ε
x
ε
y
ε
z
ε
r
x
ε
r
y
ε
r
z
]
T
,
where δ
φ,  δ
λ are the position errors of INS; δ
V
x, δ
V
y are the velocity errors of INS; α, β and γ are attitude errors of INS; ε
x, ε
y and ε
z are the gyro constant drifts; ε
r
x, ε
r
y and ε
r
z are the gyro random drifts.The state vector of GPS is selected as(14)
X
G
=
[
δ
φ
G
δ
λ
G
]
T
,
where δ
φ
G, δ
λ
G are the position error of GPS.Position integrated mode is used for INS/GPS subfilter, so the state equation and measurement equation are given by(15)
[
X
˙
I
X
˙
G
]
=
[
F
I
0
0
F
G
]
[
X
I
X
G
]
+
[
W
I
W
G
]
,
Z
I
G
=
[
φ
I
-
φ
G
λ
I
-
λ
G
]
=
H
I
G
[
X
I
X
G
]
+
V
I
G
,
where F
I is INS state transition matrix in [15]; W
I, W
G are the system noise matrixes for INS and GPS, respectively; V
I
G is measurement noise matrix; and H
I
G is state transition matrix for INS/GPS measurement equation. Consider
(16)
W
I
=
[
ω
g
x
ω
g
y
ω
g
z
]
T
,
W
G
=
[
ω
φ
ω
λ
]
T
,
where ω
g
x, ω
g
y and ω
g
z are system white noises of INS and ω
φ, ω
λ are system white noise of GPS. Consider
(17)
H
I
G
=
[
1
0
0
0
0
0
0
0
0
0
0
0
0
-
1
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
-
1
]
.Using first order Markov process to complete the approximate fitting, GPS state transition matrixF
G can be expressed as follows:
(18)
F
G
=
diag
⁡
(
1
τ
GPS
φ
,
1
τ
GPS
λ
)
,
where τ
GPS
φ and τ
GPS
λ are the correlation times.Velocity matching method is used in INS/DVL subfilter, and the state equation and measurement equation are given by(19)
[
X
˙
I
X
˙
D
]
=
[
F
I
0
0
F
D
]
[
X
I
X
D
]
+
[
W
I
W
D
]
,
Z
I
D
=
[
V
I
x
-
V
D
x
V
I
y
-
V
D
y
]
=
H
I
D
[
X
I
X
D
]
+
V
I
D
,
where X
D
=
[
δ
V
D
x
δ
V
D
y
δ
k
]
T is state vector of DVL velocity estimation error, δ
V
D
x, δ
V
D
y are velocity error, and δ
k is scale error. F
D is DVL state transition matrix:
(20)
F
D
=
diag
⁡
(
1
τ
D
x
,
1
τ
D
y
,
1
τ
k
)
,
where τ
D
x, τ
D
y and τ
k are Markov correlation time, W
D is system noise matrix, ω
D
x, ω
D
y are driven white noise, and V
I
D is measure noise matrix. Consider
(21)
W
G
=
[
ω
D
x
ω
D
y
]
T
,
H
I
D
=
[
0
0
1
0
0
0
-
V
y
0
0
0
0
0
0
-
1
0
-
V
x
0
0
0
1
0
0
V
x
0
0
0
0
0
0
0
-
1
-
V
y
]
.


## 4.2. INS/GPS/DVL Global Filter

To make the integrated filter realizable for computer calculation, it is necessary to discretize the state equation and measurement equation mentioned above and then we have(22)
X
k
=
Φ
k
,
k
-
1
X
k
+
Γ
k
-
1
W
k
-
1
,
Z
k
=
H
k
X
k
+
V
k
.Assuming that the local estimates are unrelated, the global estimate can be expressed as(23)
X
^
g
=
P
g
(
P
11
-
1
X
^
1
+
P
22
-
1
X
^
2
)
,
P
g
=
(
P
11
-
1
+
P
22
-
1
)
-
1
.Local estimates (X
^
1
,
X
^
2
) from the subfilters and their covariance matrix (P
1
,
P
2
) are sent to master filter, and the obtained results would integrate with master filter estimates to get the global optimal estimates. Besides, global estimate X
^
g  and its covariance matrix P
g, which are amplified to be β
i
-
1
P
g
(
β
i
≤
1
), are feed back to sub filters, and estimates of subfilter are reset as
(24)
X
^
i
=
X
^
g
,
P
i
i
=
β
i
-
1
P
g
.


## 5. Simulation

### 5.1. Simulation for Fault Positioning

Compared with the detection step error, slope error is more difficult to be detected. Consequently, the slope error is selected to be researched to verify the efficience of the proposed fault detection method.(1)
Simulation conditions:
During 100~200 sampling points, the slope error is added to the sensor signal.
(2)
Simulation results:As a comparison, the fault detection method based on Gaussian wavelet transform is applied to faults detection; its fault detection result is shown in Figure4.Figure 4
Slope error detection based on Gaussian wavelet transform.It is obvious that there is a time delay for the fault detection algorithm based on Gaussian wavelet transform in Figure4. And the sensor faults cannot be efficiently detected until the soft faults accumulate to a certain degree.And the simulation result of the proposed fault detection method based on binary wavelet transform and modulus maxima is shown in Figure5.Figure 5
Slope error detection by binary wavelet transform.As known from Figure5, the proposed fault detection method can efficiently detect the slope error of sensor fault and exhibits excellent accuracy for the fault positioning and sensational resolution in time domain. However, with the increase of the scale, multiple extrema appear in the signal which are caused by noise, and some of extrema would attenuate sharply or even disappear. However, those extrema which have appeared in all scales would not attenuate. And they are corresponding to the real singularities of signal.


### 5.2. Simulation for Fault Type Detection

The step fault and slope fault are selected to analyze the fault type detection performance of the proposed fault detection method.During 100~200 sampling points, the constant error (to form the step fault) and slope error are added in the sensor signals, respectively. The logarithm of modulus maxima is calculated atu
=
100, as shown in Figure 6.Figure 6
Logarithm curves of modulus maxima atu
=
100.In Figure6, the logarithm curves of modulus maxima exhibit obvious attenuation ratio at u
=
100. In order to identify the fault type, the attenuate ratio for logarithm curves of modulus maxima with respect to the scales at u
=
100 should be obtained. By using linear least square method, the slope of curve on each scale of the curve can be estimated, and the Lipschitz exponent can be obtained by subtracting 0.5 from the slope. As for the two kinds of faults, the scopes of modulus maxima curves and Lipschitz exponents at u
=
100 are shown in Table 1.Table 1
Lipschitz exponentα.
Step fault
Slope fault
Theoretical value ofα
0
1
Modulus maxima attenuation ratio
0.4737
1.5143
Measured value ofα
−0.0263
1.0143According to Table1, the measured Lipschitz exponents with respect to step fault and slope fault approximate to their ideal values, respectively. By calculating the Lipschitz exponent of signal, the fault type can be detected.In conclusion, by searching wavelet modulus maxima on different scales, the fault detection method based on modulus maxima can accurately position fault points and distinguish the fault type through Lipschitz exponent which is calculated by the attenuation of wavelet modulus maxima.


### 5.3. Simulation for Fault-Tolerant Filter

(1)
Assuming that initial latitude isφ
0
=
3
0
° and initial longitude is λ
0
=
12
0
°, gyro constant drift is 0.0
1
°
/
h; heading angel is 4
5
° and east velocity and north velocity are 5 kn; simulation time is 5.5 h.
(2)
INS/DVL/GPS integrated system is applied to calibration. After 4-hour system work, inertial system is calibrated by the gyro drift estimated and system error information.
(3)
Two conditions are involved in the calibration simulation: (a) GPS without any fault; (b) GPS with slope fault.For the two conditions, the drift estimated results for east and azimuth gyroscopes are given in Figure7.Gyro drifts estimation result.
(a)
The drift estimation results when GPS is without faults
(b)
The drift estimation results when GPS has slope faultsIt is obvious from Figure7 that on the no-faults condition the gyro drifts and position information can be well estimated, especially since the estimation error for east gyro drift is less than 10%. On the other hand, when slope fault occurs, we cannot obtain the position measurement, and the gyro drift estimation performance decreases, especially since the estimate accuracy of the azimuth gyro drift is seriously restricted by limited available information. Through the detection and isolation of slope faults in INS/GPS subintegrated system, such a single system fault will have no effect on the output of SINS/DVL system but slightly decrease for the position estimation comparing with the situation without GPS faults.


## 5.1. Simulation for Fault Positioning

Compared with the detection step error, slope error is more difficult to be detected. Consequently, the slope error is selected to be researched to verify the efficience of the proposed fault detection method.(1)
Simulation conditions:
During 100~200 sampling points, the slope error is added to the sensor signal.
(2)
Simulation results:As a comparison, the fault detection method based on Gaussian wavelet transform is applied to faults detection; its fault detection result is shown in Figure4.Figure 4
Slope error detection based on Gaussian wavelet transform.It is obvious that there is a time delay for the fault detection algorithm based on Gaussian wavelet transform in Figure4. And the sensor faults cannot be efficiently detected until the soft faults accumulate to a certain degree.And the simulation result of the proposed fault detection method based on binary wavelet transform and modulus maxima is shown in Figure5.Figure 5
Slope error detection by binary wavelet transform.As known from Figure5, the proposed fault detection method can efficiently detect the slope error of sensor fault and exhibits excellent accuracy for the fault positioning and sensational resolution in time domain. However, with the increase of the scale, multiple extrema appear in the signal which are caused by noise, and some of extrema would attenuate sharply or even disappear. However, those extrema which have appeared in all scales would not attenuate. And they are corresponding to the real singularities of signal.


## 5.2. Simulation for Fault Type Detection

The step fault and slope fault are selected to analyze the fault type detection performance of the proposed fault detection method.During 100~200 sampling points, the constant error (to form the step fault) and slope error are added in the sensor signals, respectively. The logarithm of modulus maxima is calculated atu
=
100, as shown in Figure 6.Figure 6
Logarithm curves of modulus maxima atu
=
100.In Figure6, the logarithm curves of modulus maxima exhibit obvious attenuation ratio at u
=
100. In order to identify the fault type, the attenuate ratio for logarithm curves of modulus maxima with respect to the scales at u
=
100 should be obtained. By using linear least square method, the slope of curve on each scale of the curve can be estimated, and the Lipschitz exponent can be obtained by subtracting 0.5 from the slope. As for the two kinds of faults, the scopes of modulus maxima curves and Lipschitz exponents at u
=
100 are shown in Table 1.Table 1
Lipschitz exponentα.
Step fault
Slope fault
Theoretical value ofα
0
1
Modulus maxima attenuation ratio
0.4737
1.5143
Measured value ofα
−0.0263
1.0143According to Table1, the measured Lipschitz exponents with respect to step fault and slope fault approximate to their ideal values, respectively. By calculating the Lipschitz exponent of signal, the fault type can be detected.In conclusion, by searching wavelet modulus maxima on different scales, the fault detection method based on modulus maxima can accurately position fault points and distinguish the fault type through Lipschitz exponent which is calculated by the attenuation of wavelet modulus maxima.


## 5.3. Simulation for Fault-Tolerant Filter

(1)
Assuming that initial latitude isφ
0
=
3
0
° and initial longitude is λ
0
=
12
0
°, gyro constant drift is 0.0
1
°
/
h; heading angel is 4
5
° and east velocity and north velocity are 5 kn; simulation time is 5.5 h.
(2)
INS/DVL/GPS integrated system is applied to calibration. After 4-hour system work, inertial system is calibrated by the gyro drift estimated and system error information.
(3)
Two conditions are involved in the calibration simulation: (a) GPS without any fault; (b) GPS with slope fault.For the two conditions, the drift estimated results for east and azimuth gyroscopes are given in Figure7.Gyro drifts estimation result.
(a)
The drift estimation results when GPS is without faults
(b)
The drift estimation results when GPS has slope faultsIt is obvious from Figure7 that on the no-faults condition the gyro drifts and position information can be well estimated, especially since the estimation error for east gyro drift is less than 10%. On the other hand, when slope fault occurs, we cannot obtain the position measurement, and the gyro drift estimation performance decreases, especially since the estimate accuracy of the azimuth gyro drift is seriously restricted by limited available information. Through the detection and isolation of slope faults in INS/GPS subintegrated system, such a single system fault will have no effect on the output of SINS/DVL system but slightly decrease for the position estimation comparing with the situation without GPS faults.


## 6. Conclusion

This current work has presented a sensor fault detection method based on binary wavelet transform, modulus maxima, and Lipschitz exponent. In this method, after binary wavelet transform of original signal, sensor faults can be accurately positioned by searching modulus maxima, and the type of fault can be obtained by calculation of the Lipschitz exponents. To guarantee the performance during calibration phase for INS navigation errors, a developed fault-tolerant federated filter for INS/GPS/DVL integrated navigation system is proposed. Simulation shows that even when GPS slope faults occur, INS/GPS/DVL system could keep on providing reliable position information only with a slight drop of positioning accuracy.



---

*Source: 825726-2014-08-17.xml*