# An Aeromagnetic Compensation Algorithm for Aircraft Based on Fuzzy Adaptive Kalman Filter

Authors: Yao Zhang; Yu-Xin Zhao; Shuai Chang
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

In order to ensure the effectiveness of geomagnetic navigation, as the foundation, the precise measurement of geomagnetic field must be guaranteed; namely, aircraft aeromagnetic compensation is worthy of being further studied. In this paper, the classical aircraft aeromagnetic compensation algorithm based on Leliak Model is analyzed and an aircraft aeromagnetic compensation algorithm based on fuzzy adaptive Kalman filter is proposed, which is a new approach for aircraft to achieve aeromagnetic compensation. Simulation results show that it has better compensation performance without relying on the aircraft attitude.


---

## Body

## 1. Introduction

Measurement of geomagnetic field refers to obtaining the field strength of carrier location accurately and in a real-time manner, and to improve the accuracy of measurement, there are mainly two ways: on one hand, improving the measurement precision of the sensor and, on the other hand, overcoming the interference caused by carrier’s own magnetic field [1]. At present, the research on geomagnetic navigation and measurement of geomagnetic field generally uses orthotropic three-axis magnetic sensor as the basic magnetic components, because of its unique advantages that it can measure the overall field strength under any angle conveniently and need not consider the aircraft attitude and local magnetic field vector in an ideal world [2]. However, in practical application, due to the limitation of processing and assembling accuracy level, low cost of three-axis magnetic sensor cannot ensure that three-axis is orthogonal strictly, and the electric performance of the three axis sensitivity cannot completely be symmetrical [3]. The method of eliminating the magnetic field interference is generally neglected in the current research, especially, the interference caused by carriers themselves, and this kind of interference which is relative to carrier attitude is called aeromagnetic interference [4]. Tolles and Lawson reported the existence of interference when aircraft was engine-driven in 1950, and consider it as a kind of constant noise [5]. The amplitude of aeromagnetic interference is small, which is composed by constant magnetic interference, excitation magnetic interference, and eddy current interference; however, it cannot be ignored relative to the weak magnetic field, and cannot simply it into a white noise because different aircraft attitude will cause different amplitude [6]. Therefore, it is necessary to consider the noise as a kind of complex form.Leliak had a further study on this problem and found a 16th-order linear model, called “Leliak Model”, to achieve the aeromagnetic compensation, which is used as the standard method in practical application [7]. Many scholars have studied the “Leliak Model” and improved it on certain aspect [8–11]. Bickel analyzed the magnetic interference under the small signal compensation and proposed a calculating method which is proved to be effective [8]; however, this research has bigger limitations and it cannot work when the signal is not small enough. Vijay Gopal et al. had a research on the high resolution airborne magnetic surveys and paid more attention to the real time compensation, which has better measurability and high value of application [9], but when the aircraft’s location is less than 1000 km, this method will lose precision more or less. Jia and EiKon proposed a new aircraft compensation system for magnetic terrains, which is based on the Ontario mineral exploration program, and its research object is mainly about Ontario mineral [10]; Above all, it can be seen that, by using this method to solve other field’s magnetic interference problems, it has shown limitations. Liu had a comprehensive study on aeromagnetic compensation technology including setting up interference model, proposing compensation algorithm, and providing measuring method in practical experiment [11]; however, there are some aspects to be improved, such that the measure accuracy is not high enough and the method is effective when the geomagnetic field strength can be considered as a constant only. Therefore, it is necessary to provide a new aeromagnetic compensation algorithm which has no limitation in signal’s value, aircraft flying height, and the change rate of geomagnetic field strength.Kalman filter has been widely used to forecast system state vectors and estimate system parameters since it was invented by Stanley Schmidt because of its high accuracy with strong ability of interference suppression [12]. Under the situation that the calculation errors are ignored, classical Kalman filter can obtain the state vector value accurately; however, this situation is too hard to be satisfied, and it will cause divergent phenomenon sometime [13]. To solve the problem above, use measurement information constantly online to estimate the statistical feature of measurement noise in order to get accurate estimates of system state variables when filtering. At present, there are a lot of adaptive filtering methods, and the combination of Kalman filtering and fuzzy logic method is an ideal approach [14]. Bai et al. proposed an adaptive fuzzy Kalman filter to estimate system state vector, which adjusts estimation results based on observing noise in every step [15]; however, this method will cause oscillation of the covariance matrix.Due to the fact that the existing works on aeromagnetic compensation were all found to rely on the aircraft attitude angle, the accuracy and singularity will be influenced. Therefore, it is necessary to make researches on aeromagnetic compensation algorithm using Kalman filter, which has better compensation performance without relying on the aircraft attitude.At present, the research on aeromagnetic compensation algorithm using Kalman filter is seldom. In this paper, we proposed a new way to achieve aircraft aeromagnetic compensation, which can avoid the classical method limitation of accuracy and singularity. The rest of this paper consists of four sections. Leliak Model is analyzed in Section2, and its shortcomings are pointed out. Section 3 improves the Leliak and gives the aircraft aeromagnetic compensation algorithm based on Kalman filter and adaptive fuzzy Kalman filter. Section 4 validates the effectiveness of the proposed method. Finally, we conclude this paper in Section 5.


## 2. Problem Description and Preliminaries

The classical aircraft aeromagnetic compensation algorithm based on Leliak Model can be expressed as follows:(1)
Δ
T
→
=
T
→
p
+
T
→
i
+
T
→
e
,
where T
→
p, T
→
i, and T
→
e are constant magnetic interference, excitation magnetic interference, and eddy current interference, and Δ
T
→ is the total geomagnetic interference caused by carrier moving.Constant magnetic interference is caused by permanent magnet and ferromagnetic materials magnetized permanently, which cannot change according to the carrier attitude and can be considered as a constant, and it can be divided as follows:(2)
T
p
d
=
T
p
x
cos
⁡
X
+
T
p
y
cos
⁡
Y
+
T
p
z
cos
⁡
Z
,
where T
p
d is the sensor effect caused by constant magnetic interference, T
p
x, T
p
y, and T
p
z are the projection weight of T
p
d in the carrier coordinate system, and X, Y, and Z are the angle between magnetic field direction and O
x, O
y, and O
z of carrier coordinate system.Excitation magnetic interference is caused by soft iron material magnetized by earth’s magnetic field temporarily, which is proportional with the earth’s magnetic field, and its expression is(3)
T
i
d
=
T
[
(
a
11
cos
⁡
X
+
a
12
cos
⁡
Y
+
a
13
cos
⁡
Z
)
cos
⁡
X
+
(
a
21
cos
⁡
X
+
a
22
cos
⁡
Y
+
a
23
cos
⁡
Z
)
cos
⁡
Y
+
(
a
31
cos
⁡
X
+
a
32
cos
⁡
Y
+
a
33
cos
⁡
Z
)
cos
⁡
Z
]
,
where T
i
d is the sensor effect caused by excitation magnetic interference, T is magnetic field strength, and a
i
j
(
i
,
j
=
1,2
,
3
) is the proportionality coefficient, which is only related to carrier’s material structure and material and can be considered as constants.Eddy current interference is caused by the changes of magnetic flux when the carrier’s attitude is changing, and there is eddy current among the conductor materials, which will bring interference to the magnetic field. This kind of interference is relative to the change rate of magnetic field, and its expression is(4)
T
ed
=
T
[
(
b
11
-
b
22
)
cos
⁡
X
(
-
sin
X
X
˙
)
+
b
12
cos
⁡
X
(
-
sin
Y
Y
˙
)
+
b
13
cos
⁡
X
(
-
sin
Z
Z
˙
)
]
+
b
21
cos
⁡
Y
(
-
sin
X
X
˙
)
+
b
23
cos
⁡
Y
(
-
sin
Z
Z
˙
)
+
b
32
cos
⁡
Z
(
-
sin
Y
Y
˙
)
+
(
b
33
-
b
22
)
cos
⁡
Z
(
-
sin
Z
Z
˙
)
,
where T
ed is the sensor effect caused by eddy current interference and b
i
j
(
i
,
j
=
1,2
,
3
) is the proportionality coefficient, which is only related to carrier’s material structure and material and can be considered as constants.The schematic diagram of the Leliak Model is shown in Figure1.Figure 1
The schematic diagram of the Leliak Model.From (1), (2), (3), and (4), if X, Y, Z, and T are known, we can calculate the sum interference caused by carrier moving by obtaining 16 parameters including T
p
x, T
p
y, T
p
z, (
a
11
-
a
22
), (
a
33
-
a
22
), (
a
21
+
a
12
), (
a
31
+
a
13
), (
a
32
+
a
23
), (
b
11
-
b
22
), (
b
33
-
b
22
), b
12, b
13, b
21, b
23, b
31, and b
32. Then, we can achieve the reverse compensation of the sensor measuring value according the calculation of interference value. In order to simplify the calculating process, Leliak Model uses reducing variable number method, whose specific approach is letting the aircraft do the single shaft motion to ignore the change rate of the other two attitudes.Above all, classical aircraft aeromagnetic compensation algorithm based on Leliak Model is effective to achieve the aircraft aeromagnetic compensation; however, there is coupling among the aircraft three channel motions, which cannot be ignored, and it will influence the measure precision.


## 3. Aeromagnetic Compensation Algorithm Based on Fuzzy Adaptive Kalman Filter

### 3.1. Geomagnetic Interference Model Improving

In order to take the coupling among the aircraft three channels into consideration, we can establish the mathematical model in the carrier coordinate system according to the interference magnetic field generation mechanism in the following form [16]:
(5)
T
m
=
T
+
T
p
+
A
T
+
B
T
˙
,
where T
m
=
[
T
m
x
T
m
y
T
m
z
]
T is the sensor measurement, T
=
[
T
x
T
y
T
z
]
T is the actual magnetic field strength, T
p
=
[
T
p
x
T
p
y
T
p
z
]
T is the constant magnetic interference, and A
=
{
a
i
j
}
3
×
3
(
i
,
j
=
1,2
,
3
) and B
=
{
b
i
j
}
3
×
3
(
i
,
j
=
1,2
,
3
) are the proportionality coefficient matrices.We can obtain formula (6) by transferring carrier coordinate system to the geographic coordinate system, and A
(
t
) is the direction cosine matrix,
(6)
[
T
x
T
y
T
z
]
=
A
(
t
)
[
T
x
0
T
y
0
T
z
0
]
,
where T
x
0, T
y
0, and T
z
0 are the magnetic field strength of the shafts O
x, O
y, and O
z projection in the geographic coordinate system. Combined with (6), (5) can be written in the following form:
(7)
[
T
m
x
T
m
y
T
m
z
]
=
A
(
t
)
[
T
x
o
T
y
o
T
z
o
]
+
[
T
p
x
T
p
y
T
p
z
]
+
[
a
11
a
12
a
13
a
21
a
22
a
23
a
31
a
32
a
33
]
A
(
t
)
[
T
x
o
T
y
o
T
z
o
]
+
[
b
11
b
12
b
13
b
21
b
22
b
23
b
31
b
32
b
33
]
d
A
(
t
)
d
t
[
T
x
o
T
y
o
T
z
o
]
,
where
(8)
d
A
(
t
)
d
t
=
[
0
ω
z
-
ω
y
-
ω
z
0
ω
x
ω
y
-
ω
x
0
]
A
(
t
)
,
ω
x, ω
y, and ω
z are the attitude angular velocity around the shafts O
x, O
y, and O
z in the inertial coordinate system.


### 3.2. Kalman Filter Model Establishing

From formula (7) above, it can be seen that the geomagnetic interference model can be written as
(9)
[
T
m
x
T
m
y
T
m
z
]
-
[
T
x
T
y
T
z
]
=
[
T
p
x
T
p
y
T
p
z
]
+
[
a
11
a
12
a
13
a
21
a
22
a
23
a
31
a
32
a
33
]
[
T
x
T
y
T
z
]
+
[
b
11
b
12
b
13
b
21
b
22
b
23
b
31
b
32
b
33
]
[
T
˙
x
T
˙
y
T
˙
z
]
;
its physical meaning is that the sum interference for sensor measurement consists of constant magnetic interference, excitation magnetic interference, and eddy current interference.Let(10)
x
=
[
T
p
x
T
p
y
T
p
z
[
a
11
⋯
a
33
]
1
×
9
[
b
11
⋯
b
33
]
1
×
9
]
T
,
y
=
[
T
m
x
-
T
x
T
m
y
-
T
y
T
m
z
-
T
z
]
T
,
and formula (9) can be written as
(11)
y
=
H
(
t
)
x
,
where H
(
t
) can be obtained as follows from formula (9):(12)
H
(
t
)
=
[
1
0
0
T
x
T
y
T
z
0
0
0
0
0
0
T
˙
x
T
˙
y
T
˙
z
0
0
0
0
0
0
0
1
0
0
0
0
T
x
T
y
T
z
0
0
0
0
0
0
T
˙
x
T
˙
y
T
˙
z
0
0
0
0
0
1
0
0
0
0
0
0
T
x
T
y
T
z
0
0
0
0
0
0
T
˙
x
T
˙
y
T
˙
z
]
.Kalman filter can be expressed as formula (13) with formula (11) discretization as follows:
(13)
y
(
k
)
=
H
(
k
)
x
(
k
)
,
where H
(
k
) is in the following form:(14)
H
(
k
)
=
[
E
3
×
3
T
3
×
9
(
k
)
T
3
×
9
(
k
+
1
)
]
,
E
3
×
3
=
[
1
0
0
0
1
0
0
0
1
]
,
T
3
×
9
(
k
)
=
[
T
x
(
k
)
T
y
(
k
)
T
z
(
k
)
0
0
0
0
0
0
0
0
0
T
x
(
k
)
T
y
(
k
)
T
z
(
k
)
0
0
0
0
0
0
0
0
0
T
x
(
k
)
T
y
(
k
)
T
z
(
k
)
]
,
T
3
×
9
(
k
+
1
)
=
[
T
x
(
k
+
1
)
T
y
(
k
+
1
)
T
z
(
k
+
1
)
0
0
0
0
0
0
0
0
0
T
x
(
k
+
1
)
T
y
(
k
+
1
)
T
z
(
k
+
1
)
0
0
0
0
0
0
0
0
0
T
x
(
k
+
1
)
T
y
(
k
+
1
)
T
z
(
k
+
1
)
]
.It is obvious thatH
(
t
) can be treated as the observe matrix of the Kalman filter, the state vector is x
(
k
), and the observe vector is y
(
k
).Since the elements of state vector are constants, the state matrix of the Kalman filter is unit matrix; namely, the state equation is(15)
x
(
k
)
=
E
x
(
k
-
1
)
+
W
(
k
)
,
where E is an unit matrix, u
(
k
) is controller, and W
(
k
) is the process noise, which is white Gaussian noise and its covariance is Q.Take observe error into consideration, the observe equation is(16)
y
(
k
)
=
H
(
k
)
x
(
k
)
+
V
(
k
)
,
where V
(
k
) is the observe noise, which is white Gaussian noise and its covariance is R.Above all, Kalman filter for aeromagnetic compensation has been established which is expressed by state equation (15) and observe equation (16).


### 3.3. Aeromagnetic Compensation Algorithm Based on Kalman Filter

Aeromagnetic compensation algorithm proposed here uses Kalman filter method to solve the problem that classical aeromagnetic compensation algorithm cannot ensure the geomagnetic interference being measured accurately [17, 18], and the aeromagnetic compensation algorithm based on Kalman filter is introduced in this section. For a linear stochastic differential system, if process noise and observe noise are both white Gaussian noise, Kalman filter is the optimal method to estimate state vector. Aeromagnetic compensation algorithm based on Kalman filter is shown as Figure 2.Figure 2
Aeromagnetic compensation algorithm based on Kalman filter.Firstly, the state equation of Kalman filter is used to achieve the state predicting, and it is in the following form:(17)
x
(
k
∣
k
-
1
)
=
Φ
x
(
k
-
1
∣
k
-
1
)
+
B
u
(
k
)
,
where x
(
k
∣
k
-
1
) is the predicting result at the previous state, x
(
k
-
1
∣
k
-
1
) is the optimal result at the previous state, Φ is state transition matrix, which is treated as unit matrix, and u
(
k
) is the controller at current state, which is treated as zero in aeromagnetic compensation algorithm combined with formula (15). Namely, formula (17) can be written as x
(
k
∣
k
-
1
)
=
x
(
k
-
1
∣
k
-
1
). The current state has been updated by formula (17), but the covariance of x
(
k
∣
k
-
1
) has not been updated.Secondly, we calculate the covariance as follows:(18)
P
(
k
∣
k
-
1
)
=
Φ
P
(
k
-
1
∣
k
-
1
)
Φ
T
+
Q
,
where P
(
k
∣
k
-
1
) is the covariance of x
(
k
∣
k
-
1
) and P
(
k
-
1
∣
k
-
1
) is the covariance of x
(
k
-
1
∣
k
-
1
). As Φ
=
E
21
×
21, formula (18) can be written as P
(
k
∣
k
-
1
)
=
P
(
k
-
1
∣
k
-
1
)
+
Q. Formula (17) and formula (18) bring out the predicting result of system state.Thirdly, combined with the predicting results and measuring results, the optimal value ofx
(
k
) can be obtained as
(19)
x
(
k
∣
k
)
=
x
(
k
∣
k
-
1
)
+
K
g
(
k
)
(
y
(
k
)
-
H
(
k
)
x
(
k
∣
k
-
1
)
)
,
where K
g
=
P
(
k
∣
k
-
1
)
H
T
(
k
)
/
(
H
(
k
)
P
(
k
∣
k
-
1
)
H
T
(
k
)
+
R
) is the Kalman gain.Fourthly, in order to ensure that Kalman filter keeps working until system process is finished, the covariance ofx
(
k
∣
k
) must be updated by calculating the following formula:
(20)
P
(
k
k
)
=
(
E
-
K
g
(
k
)
H
(
k
)
)
P
(
k
k
-
1
)
.By repeating iteration process above, we can get the estimatedy
(
k
) when the filter converges under the situation that the initial x
(
k
) is small enough.Above all, the aeromagnetic compensation algorithm estimates the elements ofy
(
k
) and there is no need to ask aircraft to do single shaft motion with no errors during this process. After obtaining the elements of y
(
k
), we can bring them into aircraft magnetic field model to compensate. However, there will exist divergent phenomenon if the noise cannot be estimated accurately or the aircraft geomagnetic interference is too big to measure.


### 3.4. Aeromagnetic Compensation Algorithm Based on Fuzzy Adaptive Kalman Filter

In order to ensure the effectiveness of kalman filter under the noise unmeasured, fuzzy adaptive method is considered to estimate the process noise and observe noise, which can avoid divergent phenomenon caused by noise uncertainty [19–23].Using adaptive method to adjust process noiseQ and observe noise R, we can get the following formulas:
(21)
Q
(
k
k
)
=
s
a
(
b
a
)
Q
(
k
-
1
k
-
1
)
,
R
(
k
∣
k
)
=
s
r
(
b
r
)
R
(
k
-
1
∣
k
-
1
)
,
where s
a
(
b
a
) is the adjustment parameter of Q, s
r
(
b
r
) is the adjustment parameter of R, and b
a
>
0
,
b
r
>
0 are constants standing for zoom level of Q and R, which are calculated by fuzzy inference system (FIS). Then, analysis of b
a and b
r is given as follows [24].(i) Considerb
a
,
b
r
>
1; magnify the effect on Q and R by s
a
(
b
a
) and s
r
(
b
r
). In this case, Q and R can approach the real process noise and observe noise in fewer steps. If the values of b
a
,
b
r are too big, the approaching rate is fast enough; however, it may cause small amplitude oscillation around the real value of noise.(ii) Consider0
<
b
a
,
b
r
<
1; shrink the effect on Q and R by s
a
(
b
a
) and s
r
(
b
r
). In this case, there can be no oscillation and stability at the real noise value; however, it will cause more approach steps; namely, the adjustment time will be too long.(iii) Considerb
a
,
b
r
=
0; there is no effect on Q and R by s
a
(
b
a
) and s
r
(
b
r
). In this case, this method is equal to the method proposed in Section 3.3.In order to establish the FIS of aeromagnetic compensation algorithm, performance indicators used to evaluate estimation precision are introduced here with its definition as(22)
R
O
R
(
k
)
=
T
r
(
c
r
)
T
r
(
p
r
)
,
Q
O
Q
(
k
)
=
T
r
(
c
a
)
T
r
(
p
a
)
,
where T
r
(
·
) is the trace of the matrix, p
r and p
a are residual variance of process noise and observe noise, and c
r and c
a are variance of the real value of process noise and observe noise, and their definitions are as follows:
(23)
p
r
=
H
(
k
)
(
P
(
k
-
1
∣
k
-
1
)
+
Q
(
k
∣
k
)
)
H
T
(
k
)
+
R
(
k
-
1
∣
k
-
1
)
,
p
q
=
P
(
k
-
1
∣
k
-
1
)
+
R
(
k
∣
k
)
+
Q
(
k
-
1
∣
k
-
1
)
,
(24)
c
r
=
c
q
=
1
M
∑
i
=
i
0
k
r
i
r
i
T
,
(25)
r
k
=
y
(
k
)
-
H
(
k
)
x
(
k
∣
k
-
1
)
,
where r
k is residual error sequence.From formulas (22) and (23), it can be seen that if the noise model is accurate enough, R
O
R
(
k
)
=
1
,
Q
O
Q
(
k
)
=
1 can be established, and we use R
(
k
-
1
∣
k
-
1
) to ensure R
O
R
(
k
) near to 1 and use Q
(
k
-
1
∣
k
-
1
) to ensure Q
O
Q
(
k
) near to 1.Then, using FIS to calculates
a
(
b
a
) and s
r
(
b
r
) with the following FIS rules.(i) For observe noise, the FIS rules are as follows.IfR
O
R
(
k
)
∈
{
|
R
O
R
(
k
)
-
1
|
≤
ε
}, then s
r
(
b
r
)
∈
{
|
s
r
(
b
r
)
-
1
|
≤
ε
}.
IfR
O
R
(
k
)
∈
{
R
O
R
(
k
)
>
1
+
ε
}, then s
r
(
b
r
)
∈
{
s
r
(
b
r
)
>
1
+
ε
}.
IfR
O
R
(
k
)
∈
{
R
O
R
(
k
)
<
1
-
ε
}, then s
r
(
b
r
)
∈
{
s
r
(
b
r
)
<
1
-
ε
},where0
<
ε
<
1 and ε stands for boundary range of observe noise FIS.(ii) For process noise, the FIS rules are as follows.IfQ
O
Q
(
k
)
∈
{
|
Q
O
Q
(
k
)
-
1
|
≤
η
}, then s
a
(
b
a
)
∈
{
|
s
a
(
b
a
)
-
1
|
≤
η
}.
IfQ
O
Q
(
k
)
∈
{
Q
O
Q
(
k
)
>
1
+
η
}, then s
a
(
b
a
)
∈
{
s
a
(
b
a
)
>
1
+
η
}.
IfQ
O
Q
(
k
)
∈
{
Q
O
Q
(
k
)
<
1
-
η
}, then s
a
(
b
a
)
∈
{
s
a
(
b
a
)
<
1
-
η
},where0
<
η
<
1 and η stands for boundary range of process noise FIS.The output of FIS can be plugged in formula (21) to obtain the present noise value of R
(
k
∣
k
) and Q
(
k
∣
k
), and we achieve the process noise and observe noise estimation by FIS, which can avoid divergent phenomenon in theory.


## 3.1. Geomagnetic Interference Model Improving

In order to take the coupling among the aircraft three channels into consideration, we can establish the mathematical model in the carrier coordinate system according to the interference magnetic field generation mechanism in the following form [16]:
(5)
T
m
=
T
+
T
p
+
A
T
+
B
T
˙
,
where T
m
=
[
T
m
x
T
m
y
T
m
z
]
T is the sensor measurement, T
=
[
T
x
T
y
T
z
]
T is the actual magnetic field strength, T
p
=
[
T
p
x
T
p
y
T
p
z
]
T is the constant magnetic interference, and A
=
{
a
i
j
}
3
×
3
(
i
,
j
=
1,2
,
3
) and B
=
{
b
i
j
}
3
×
3
(
i
,
j
=
1,2
,
3
) are the proportionality coefficient matrices.We can obtain formula (6) by transferring carrier coordinate system to the geographic coordinate system, and A
(
t
) is the direction cosine matrix,
(6)
[
T
x
T
y
T
z
]
=
A
(
t
)
[
T
x
0
T
y
0
T
z
0
]
,
where T
x
0, T
y
0, and T
z
0 are the magnetic field strength of the shafts O
x, O
y, and O
z projection in the geographic coordinate system. Combined with (6), (5) can be written in the following form:
(7)
[
T
m
x
T
m
y
T
m
z
]
=
A
(
t
)
[
T
x
o
T
y
o
T
z
o
]
+
[
T
p
x
T
p
y
T
p
z
]
+
[
a
11
a
12
a
13
a
21
a
22
a
23
a
31
a
32
a
33
]
A
(
t
)
[
T
x
o
T
y
o
T
z
o
]
+
[
b
11
b
12
b
13
b
21
b
22
b
23
b
31
b
32
b
33
]
d
A
(
t
)
d
t
[
T
x
o
T
y
o
T
z
o
]
,
where
(8)
d
A
(
t
)
d
t
=
[
0
ω
z
-
ω
y
-
ω
z
0
ω
x
ω
y
-
ω
x
0
]
A
(
t
)
,
ω
x, ω
y, and ω
z are the attitude angular velocity around the shafts O
x, O
y, and O
z in the inertial coordinate system.


## 3.2. Kalman Filter Model Establishing

From formula (7) above, it can be seen that the geomagnetic interference model can be written as
(9)
[
T
m
x
T
m
y
T
m
z
]
-
[
T
x
T
y
T
z
]
=
[
T
p
x
T
p
y
T
p
z
]
+
[
a
11
a
12
a
13
a
21
a
22
a
23
a
31
a
32
a
33
]
[
T
x
T
y
T
z
]
+
[
b
11
b
12
b
13
b
21
b
22
b
23
b
31
b
32
b
33
]
[
T
˙
x
T
˙
y
T
˙
z
]
;
its physical meaning is that the sum interference for sensor measurement consists of constant magnetic interference, excitation magnetic interference, and eddy current interference.Let(10)
x
=
[
T
p
x
T
p
y
T
p
z
[
a
11
⋯
a
33
]
1
×
9
[
b
11
⋯
b
33
]
1
×
9
]
T
,
y
=
[
T
m
x
-
T
x
T
m
y
-
T
y
T
m
z
-
T
z
]
T
,
and formula (9) can be written as
(11)
y
=
H
(
t
)
x
,
where H
(
t
) can be obtained as follows from formula (9):(12)
H
(
t
)
=
[
1
0
0
T
x
T
y
T
z
0
0
0
0
0
0
T
˙
x
T
˙
y
T
˙
z
0
0
0
0
0
0
0
1
0
0
0
0
T
x
T
y
T
z
0
0
0
0
0
0
T
˙
x
T
˙
y
T
˙
z
0
0
0
0
0
1
0
0
0
0
0
0
T
x
T
y
T
z
0
0
0
0
0
0
T
˙
x
T
˙
y
T
˙
z
]
.Kalman filter can be expressed as formula (13) with formula (11) discretization as follows:
(13)
y
(
k
)
=
H
(
k
)
x
(
k
)
,
where H
(
k
) is in the following form:(14)
H
(
k
)
=
[
E
3
×
3
T
3
×
9
(
k
)
T
3
×
9
(
k
+
1
)
]
,
E
3
×
3
=
[
1
0
0
0
1
0
0
0
1
]
,
T
3
×
9
(
k
)
=
[
T
x
(
k
)
T
y
(
k
)
T
z
(
k
)
0
0
0
0
0
0
0
0
0
T
x
(
k
)
T
y
(
k
)
T
z
(
k
)
0
0
0
0
0
0
0
0
0
T
x
(
k
)
T
y
(
k
)
T
z
(
k
)
]
,
T
3
×
9
(
k
+
1
)
=
[
T
x
(
k
+
1
)
T
y
(
k
+
1
)
T
z
(
k
+
1
)
0
0
0
0
0
0
0
0
0
T
x
(
k
+
1
)
T
y
(
k
+
1
)
T
z
(
k
+
1
)
0
0
0
0
0
0
0
0
0
T
x
(
k
+
1
)
T
y
(
k
+
1
)
T
z
(
k
+
1
)
]
.It is obvious thatH
(
t
) can be treated as the observe matrix of the Kalman filter, the state vector is x
(
k
), and the observe vector is y
(
k
).Since the elements of state vector are constants, the state matrix of the Kalman filter is unit matrix; namely, the state equation is(15)
x
(
k
)
=
E
x
(
k
-
1
)
+
W
(
k
)
,
where E is an unit matrix, u
(
k
) is controller, and W
(
k
) is the process noise, which is white Gaussian noise and its covariance is Q.Take observe error into consideration, the observe equation is(16)
y
(
k
)
=
H
(
k
)
x
(
k
)
+
V
(
k
)
,
where V
(
k
) is the observe noise, which is white Gaussian noise and its covariance is R.Above all, Kalman filter for aeromagnetic compensation has been established which is expressed by state equation (15) and observe equation (16).


## 3.3. Aeromagnetic Compensation Algorithm Based on Kalman Filter

Aeromagnetic compensation algorithm proposed here uses Kalman filter method to solve the problem that classical aeromagnetic compensation algorithm cannot ensure the geomagnetic interference being measured accurately [17, 18], and the aeromagnetic compensation algorithm based on Kalman filter is introduced in this section. For a linear stochastic differential system, if process noise and observe noise are both white Gaussian noise, Kalman filter is the optimal method to estimate state vector. Aeromagnetic compensation algorithm based on Kalman filter is shown as Figure 2.Figure 2
Aeromagnetic compensation algorithm based on Kalman filter.Firstly, the state equation of Kalman filter is used to achieve the state predicting, and it is in the following form:(17)
x
(
k
∣
k
-
1
)
=
Φ
x
(
k
-
1
∣
k
-
1
)
+
B
u
(
k
)
,
where x
(
k
∣
k
-
1
) is the predicting result at the previous state, x
(
k
-
1
∣
k
-
1
) is the optimal result at the previous state, Φ is state transition matrix, which is treated as unit matrix, and u
(
k
) is the controller at current state, which is treated as zero in aeromagnetic compensation algorithm combined with formula (15). Namely, formula (17) can be written as x
(
k
∣
k
-
1
)
=
x
(
k
-
1
∣
k
-
1
). The current state has been updated by formula (17), but the covariance of x
(
k
∣
k
-
1
) has not been updated.Secondly, we calculate the covariance as follows:(18)
P
(
k
∣
k
-
1
)
=
Φ
P
(
k
-
1
∣
k
-
1
)
Φ
T
+
Q
,
where P
(
k
∣
k
-
1
) is the covariance of x
(
k
∣
k
-
1
) and P
(
k
-
1
∣
k
-
1
) is the covariance of x
(
k
-
1
∣
k
-
1
). As Φ
=
E
21
×
21, formula (18) can be written as P
(
k
∣
k
-
1
)
=
P
(
k
-
1
∣
k
-
1
)
+
Q. Formula (17) and formula (18) bring out the predicting result of system state.Thirdly, combined with the predicting results and measuring results, the optimal value ofx
(
k
) can be obtained as
(19)
x
(
k
∣
k
)
=
x
(
k
∣
k
-
1
)
+
K
g
(
k
)
(
y
(
k
)
-
H
(
k
)
x
(
k
∣
k
-
1
)
)
,
where K
g
=
P
(
k
∣
k
-
1
)
H
T
(
k
)
/
(
H
(
k
)
P
(
k
∣
k
-
1
)
H
T
(
k
)
+
R
) is the Kalman gain.Fourthly, in order to ensure that Kalman filter keeps working until system process is finished, the covariance ofx
(
k
∣
k
) must be updated by calculating the following formula:
(20)
P
(
k
k
)
=
(
E
-
K
g
(
k
)
H
(
k
)
)
P
(
k
k
-
1
)
.By repeating iteration process above, we can get the estimatedy
(
k
) when the filter converges under the situation that the initial x
(
k
) is small enough.Above all, the aeromagnetic compensation algorithm estimates the elements ofy
(
k
) and there is no need to ask aircraft to do single shaft motion with no errors during this process. After obtaining the elements of y
(
k
), we can bring them into aircraft magnetic field model to compensate. However, there will exist divergent phenomenon if the noise cannot be estimated accurately or the aircraft geomagnetic interference is too big to measure.


## 3.4. Aeromagnetic Compensation Algorithm Based on Fuzzy Adaptive Kalman Filter

In order to ensure the effectiveness of kalman filter under the noise unmeasured, fuzzy adaptive method is considered to estimate the process noise and observe noise, which can avoid divergent phenomenon caused by noise uncertainty [19–23].Using adaptive method to adjust process noiseQ and observe noise R, we can get the following formulas:
(21)
Q
(
k
k
)
=
s
a
(
b
a
)
Q
(
k
-
1
k
-
1
)
,
R
(
k
∣
k
)
=
s
r
(
b
r
)
R
(
k
-
1
∣
k
-
1
)
,
where s
a
(
b
a
) is the adjustment parameter of Q, s
r
(
b
r
) is the adjustment parameter of R, and b
a
>
0
,
b
r
>
0 are constants standing for zoom level of Q and R, which are calculated by fuzzy inference system (FIS). Then, analysis of b
a and b
r is given as follows [24].(i) Considerb
a
,
b
r
>
1; magnify the effect on Q and R by s
a
(
b
a
) and s
r
(
b
r
). In this case, Q and R can approach the real process noise and observe noise in fewer steps. If the values of b
a
,
b
r are too big, the approaching rate is fast enough; however, it may cause small amplitude oscillation around the real value of noise.(ii) Consider0
<
b
a
,
b
r
<
1; shrink the effect on Q and R by s
a
(
b
a
) and s
r
(
b
r
). In this case, there can be no oscillation and stability at the real noise value; however, it will cause more approach steps; namely, the adjustment time will be too long.(iii) Considerb
a
,
b
r
=
0; there is no effect on Q and R by s
a
(
b
a
) and s
r
(
b
r
). In this case, this method is equal to the method proposed in Section 3.3.In order to establish the FIS of aeromagnetic compensation algorithm, performance indicators used to evaluate estimation precision are introduced here with its definition as(22)
R
O
R
(
k
)
=
T
r
(
c
r
)
T
r
(
p
r
)
,
Q
O
Q
(
k
)
=
T
r
(
c
a
)
T
r
(
p
a
)
,
where T
r
(
·
) is the trace of the matrix, p
r and p
a are residual variance of process noise and observe noise, and c
r and c
a are variance of the real value of process noise and observe noise, and their definitions are as follows:
(23)
p
r
=
H
(
k
)
(
P
(
k
-
1
∣
k
-
1
)
+
Q
(
k
∣
k
)
)
H
T
(
k
)
+
R
(
k
-
1
∣
k
-
1
)
,
p
q
=
P
(
k
-
1
∣
k
-
1
)
+
R
(
k
∣
k
)
+
Q
(
k
-
1
∣
k
-
1
)
,
(24)
c
r
=
c
q
=
1
M
∑
i
=
i
0
k
r
i
r
i
T
,
(25)
r
k
=
y
(
k
)
-
H
(
k
)
x
(
k
∣
k
-
1
)
,
where r
k is residual error sequence.From formulas (22) and (23), it can be seen that if the noise model is accurate enough, R
O
R
(
k
)
=
1
,
Q
O
Q
(
k
)
=
1 can be established, and we use R
(
k
-
1
∣
k
-
1
) to ensure R
O
R
(
k
) near to 1 and use Q
(
k
-
1
∣
k
-
1
) to ensure Q
O
Q
(
k
) near to 1.Then, using FIS to calculates
a
(
b
a
) and s
r
(
b
r
) with the following FIS rules.(i) For observe noise, the FIS rules are as follows.IfR
O
R
(
k
)
∈
{
|
R
O
R
(
k
)
-
1
|
≤
ε
}, then s
r
(
b
r
)
∈
{
|
s
r
(
b
r
)
-
1
|
≤
ε
}.
IfR
O
R
(
k
)
∈
{
R
O
R
(
k
)
>
1
+
ε
}, then s
r
(
b
r
)
∈
{
s
r
(
b
r
)
>
1
+
ε
}.
IfR
O
R
(
k
)
∈
{
R
O
R
(
k
)
<
1
-
ε
}, then s
r
(
b
r
)
∈
{
s
r
(
b
r
)
<
1
-
ε
},where0
<
ε
<
1 and ε stands for boundary range of observe noise FIS.(ii) For process noise, the FIS rules are as follows.IfQ
O
Q
(
k
)
∈
{
|
Q
O
Q
(
k
)
-
1
|
≤
η
}, then s
a
(
b
a
)
∈
{
|
s
a
(
b
a
)
-
1
|
≤
η
}.
IfQ
O
Q
(
k
)
∈
{
Q
O
Q
(
k
)
>
1
+
η
}, then s
a
(
b
a
)
∈
{
s
a
(
b
a
)
>
1
+
η
}.
IfQ
O
Q
(
k
)
∈
{
Q
O
Q
(
k
)
<
1
-
η
}, then s
a
(
b
a
)
∈
{
s
a
(
b
a
)
<
1
-
η
},where0
<
η
<
1 and η stands for boundary range of process noise FIS.The output of FIS can be plugged in formula (21) to obtain the present noise value of R
(
k
∣
k
) and Q
(
k
∣
k
), and we achieve the process noise and observe noise estimation by FIS, which can avoid divergent phenomenon in theory.


## 4. Simulation and Analysis

In this section, experimental simulations will be carried out to evaluate the effectiveness of the proposed aeromagnetic compensation algorithm based on adaptive fuzzy Kalman filter.Detailed parameters of this algorithm are presented as follows. We put the actual magnetic field data into adaptive Kalman filter proposed in this paper, and the observe noise and process noise are both white noise with 20 nT standard deviation in each axis. Initial condition isQ
(
0
)
=
1
0
-
3
×
E
21
×
21, R
(
0
)
=
2
×
1
0
-
3
×
E
21
×
21, and x
(
0
)
=
[
0.1
0.1
0.1
[
10
⋯
10
]
1
×
9
[
15
⋯
15
]
1
×
9
]
T. System parameters are ε
=
0.05 and η
=
0.025. Simulation diagrams are as shown in Figures 3, 4, and 5.Figure 3
Magnetic field calculated by proposed aeromagnetic compensation algorithm.Kalman filter output of system state vector elements.
(a)
(b)
(c)Compensation results of three-axis magnetic field andROR,QOQ.
(a)
(b)
(c)
(d)From Figure3, we can see that the calculated magnetic field strength conforms to the distribution law of actual magnetic field, and the 60 nT total noise is measured by Kalman filter. Figure 4 shows the Kalman filter output of system state vector elements, and it can be seen that these elements convergence to constants in a finite time, which stand for the compensation result calculated by Kalman filter. From Figure 5 above, it can be seen that the noise proposed algorithm estimated is nearly equal to the actual noise in each axis; namely, the proposed aircraft aeromagnetic compensation algorithm can calculate magnetic field effectively. Furthermore, the performance indicators used to evaluate estimation precision R
O
R
(
k
) and Q
O
Q
(
k
) are around 1, which is to say that the adaptive FIS is effective in adjusting s
a
(
b
a
) and s
r
(
b
r
); namely, the variance of process noise and observe noise is estimated accurately.


## 5. Conclusions

The limitations of classical aircraft aeromagnetic compensation algorithm based on Leliak Model have been analyzed in this paper and an aircraft aeromagnetic compensation algorithm based on fuzzy adaptive Kalman filter was proposed, which provided a new approach for aircraft to achieve aeromagnetic compensation. Simulation results confirmed the compensation performance of proposed algorithm and proved that FIS is effective in estimating process noise and observe noise. Furthermore, this algorithm does not rely on the aircraft attitude, which has high application value.



---

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