# Research on the Calibration Method of Micro Inertial Measurement Unit for Engineering Application

Authors: Yun Xu; Yu Wang; Yan Su; Xinhua Zhu
Journal: Journal of Sensors
 (2016)

---

## Abstract

With rapid development of micro fabrication technology, the production level of micro inertial devices has increased, which promoted the performance improvement of Micro Inertial Measurement Unit (MIMU). Measurement precision of MIMU is one of the most significant indexes, especially for the application of the guided spinning projectiles. In order to improve the measurement precision of MIMU, this paper presents a novel calibration method. The calibration model is established and the derivation for parameters estimation has been introduced. By the multirate tests and multiposition tests, all the parameters in the calibration model can be well estimated. Verification experiment shows that the proposed method has the same compensation effect as the traditional method, but it can alleviate the computing burden for the system. Thus the proposed method will have a wide application prospect for the future engineering calibration.


---

## Body

## 1. Introduction

Nowadays, along with the development of micro fabrication technology, the manufacturing level of micro inertial devices has increased, which sprung up the production of Micro Inertial Measurement Unit (MIMU).In general, MIMU contains three Microelectromechanical Systems (MEMS) gyroscopes and accelerometers [1]. However, profit from the improvement of the integrated circuit manufacture, commercial MIMU such as HG1930, STIM300, and ADIS16488A in the market, has developed towards multiple degrees of freedom [2–4]. Digital magnetometers, pressure sensor, or some other sensors are integrated on them. Therefore, the whole performance has greatly enhanced and the aforementioned MIMU have the advantages of light weight, small size, low cost, high reliability, and long lifespan [5, 6].Generally, MEMS gyroscopes and accelerometers of MIMU are installed on a 3D orthogonal plane under the ideal condition. However, for the reason of the machining error, 3D orthogonal planes will not be completely orthogonal with each other, which will induce the installation errors for the measurement. Moreover, the bias errors and scale factor errors of the inertial sensors will also affect the measurement precision. In order to satisfy the application requirements of the designed MIMU, all of these errors should be reduced or compensated. Thus, the calibration for MIMU should be carried out.Before implementing the calibration, the calibration model should be established. As shown in the literature, the calibration model is usually established based on the bias errors, the scale factor errors, the installation errors, and some other related parameters [5, 7–9]. By means of the calibration, all parameters in the calibration model can be obtained. The calibration method for MIMU can be divided into two kinds. One is the traditional calibration method and the other is the systematic calibration method. Traditional calibration method has been well developed. It needs the high precision equipment (such as a turntable) to provide the orientation and the rate excitations, and then the parameters in the calibration model can be obtained by the optimal estimation algorithm [10–12]. Zhang et al. proposed an improved multiposition calibration method on the fact that the norms of the measured outputs of the gyroscope and accelerometer are equal to the input rotational velocity and the specific force. And then the real tests certificated that the proposed method outperformed the traditional method [13]. Traditional calibration method has the advantage of high accuracy, but the disadvantage is that the calibration precision of the parameters in the calibration model can be easily affected by the error of the turntable control system and the installation precision between MIMU and the turntable. Systematic calibration method will observe the errors of the navigation parameters (such as the attitude errors, the velocity errors, and the position errors), and then the parameters in the calibration model can be estimated by some estimation algorithms [14–17], such as fitting methods and filtering method. Fitting method should establish the particular rotation excitation for the parameters. Savge proposed two kinds of rotation procedures and then realized the estimation of the parameters in the calibration model by fitting method [18]. Filtering method often applies Kalman filter to estimate the parameters in the calibration model. By observing the velocity errors or the position errors of the navigation system under different motion, the state vector which contains the parameters in the calibration model can be estimated. Based on the principles that the magnitude of the specific force should be equal to the gravity vector and the gravity vector should be equal to the computed gravity vector which is calculated by the gyroscope outputs, Ding et al. proposed multiposition calibration method for low cost MEMS IMU and applied a Kalman filter to estimate the parameters in the calibration model [19]. Systematic calibration method has the advantage that the calibration can be implemented without any special equipment, but the disadvantage is that the calibration accuracy is limited.In the abovementioned methods, since the parameters have been obtained, all parameters should be substituted into the calibration model to get the actual input angular rates and specific forces along three axes of MIMU. During this process, the calculation of the inverse matrix [20] will be needed, which obviously increases the computing burden of the system.As we know, the bias errors, the scale factor errors, and the installation errors of the inertial sensors will affect the measurement precision of the MIMU. In this paper, we provide a new calibration method for MIMU from the point view of the engineering application. Different from the exiting method, this method does not need to calculate the inverse matrix; thus it can alleviate the computing burden of the system.For the reason that MIMU is designed to be used for the application of the guided spinning projectile, high range gyroscope is selected for the measurement of rotational speed along the projectile axis in this design. MSG7100D-3600 has the dynamic range of ±3600°/s for angular velocity, and it is installed alongz-axis to measure the rotation speed along the projectile axis. MSG7100D-300 has the dynamic range of ±300°/s, and two MSG7100D-300s are installed along x-axis and y-axis, respectively. When the projectile is lunched, there exists high shock. The shock survivability of MSA6000D-50 is 20000 g; thus three accelerometers are installed along three axes, which can bear any shock from projectile axis. The designed MIMU is illustrated in Figure 1.Figure 1
The designed MIMU.The structure of this work is as follows. In Section1, the common calibration methods for MIMU have been introduced, in which the advantages and disadvantages of them have been summarized. And then the designed MIMU which needs to be calibrated has been presented. In Section 2, the calibration method for the designed MIMU will be proposed, which contains the established calibration model, the calibration procedure, and the derivation for parameter estimation in the calibration model. In Section 3, the experimental tests will be carried out according to the given procedures in Section 2. And then error analysis of the calibration results will be discussed. Finally, conclusion will be drawn in Section 4.


## 2. Calibration Method for MIMU

### 2.1. Model Establishment for MIMU

In order to have better application for engineering application, we establish the following calibration model for the outputs of MIMU which contains the triaxial angular rates and specific forces. Different from the existing calibration model, the following model can avoid the calculation of the inverse matrix; thus it can alleviate the computing burden for the system.Output model for triaxial angular rates is as follows:(1)
ω
x
ω
y
ω
z
=
k
g
11
k
g
12
k
g
13
k
g
21
k
g
22
k
g
23
k
g
31
k
g
32
k
g
33
G
x
G
y
G
z
+
ω
x
0
ω
y
0
ω
z
0
.Output model for triaxial specific forces is as follows:(2)
f
x
f
y
f
z
=
k
a
11
k
a
12
k
a
13
k
a
21
k
a
22
k
a
23
k
a
31
k
a
32
k
a
33
A
x
A
y
A
z
+
f
x
0
f
y
0
f
z
0
,where ω
i and f
i
(
i
=
x
,
y
,
z
) mean the output angular rate and specific force of x-axis, y-axis, and z-axis, respectively. G
i and A
i
(
i
=
x
,
y
,
z
) mean the measurement of gyroscope and accelerometer; ω
i
0 and f
i
0
(
i
=
x
,
y
,
z
) mean the bias errors of the angular rate and specific force of x-axis, y-axis, and z-axis. k
i
m
n (i
=
g or a, m
,
n
=
1,2
,
3) is the coefficient related to the measurements of three gyroscopes and accelerometers, which needs to be calibrated.


### 2.2. Calibration Procedure for MIMU

Three-axis turntable applied for the calibration is shown in Figure2. Install MIMU on the plane of the inner frame, while x-axis of MIMU is along with the axis of the medium frame in the turntable, y-axis points to the front and is perpendicular to x-axis, and z-axis is along with the axis of the inner frame in the turntable. Turn the power on, and preheat MIMU for 10 minutes. And then implement the calibration procedure. Calibration procedure includes multirate tests and multiposition tests.Figure 2
Installation of MIMU on the turntable.


#### 2.2.1. Multirate Tests

After the initialization of the turntable, multirate tests will be carried out as follows.Step 1.
Adjust the medium frame in order to setx-axis point up, and acquire gyroscope outputs from MIMU when the turntable keeps stationary; denote by G
x
(
0
)
=
[
G
x
x
(
0
)
;
G
x
y
(
0
)
;
G
x
z
(
0
)
]. Rotate the outer frame with the following rotation rates: ±6°/s, ±12°/s, ±24°/s, ±30°/s, ±60°/s, ±120°/s, ±180°/s, ±240°/s, and ±300°/s, and then acquire gyroscope outputs from MIMU. If the rotation rates are positive, denote by [
G
x
x
+
(
i
)
;
G
x
y
+
(
i
)
;
G
x
z
+
(
i
)
]
(
i
=
1
~
9
). If the rotation rates are negative, denote by [
G
x
x
-
(
i
)
;
G
x
y
-
(
i
)
;
G
x
z
-
(
i
)
]
(
i
=
1
~
9
).Step 2.
Adjust the inner frame to sety-axis point up, and acquire the gyroscope outputs from MIMU when the turntable keeps stationary; denote by G
y
(
0
)
=
[
G
y
x
0
;
G
y
y
(
0
)
;
G
y
z
(
0
)
]. By the same way, rotate the outer frame with the rotation rates as in Step 1 and acquire gyroscope outputs from MIMU. If the rotation rates are positive, denote by [
G
y
x
+
(
i
)
;
G
y
y
+
(
i
)
;
G
y
z
+
(
i
)
]
(
i
=
1
~
9
). If the rotation rates are negative, denote by [
G
y
x
-
(
i
)
;
G
y
y
-
(
i
)
;
G
y
z
-
(
i
)
]
(
i
=
1
~
9
).Step 3.
Adjust the medium frame to setz-axis points up, and acquire gyroscope outputs from MIMU when the turntable keeps stationary; denote by G
z
(
0
)
=
[
G
z
x
(
0
)
;
G
z
y
(
0
)
;
G
z
z
(
0
)
]. Rotate the inner frame with the following rotation rates: 0°/s, ±10°/s, ±40°/s, ±100°/s, ±240°/s, ±400°/s, ±720°/s, ±1000°/s, ±1800°/s, ±2160°/s, ±2880°/s, and ±3600°/s, and then acquire gyroscope outputs from MIMU. If the rotation rates are positive, denote by [
G
z
x
+
(
i
)
;
G
z
y
+
(
i
)
;
G
z
z
+
(
i
)
]
(
i
=
1
~
11
). If the rotation rates are negative, denote by [
G
z
x
-
(
i
)
;
G
z
y
-
(
i
)
;
G
z
z
-
(
i
)
]
(
i
=
1
~
11
).


#### 2.2.2. Multiposition Tests

After multirate tests have been accomplished, adjust three axes of the turntable as the initialization.Step 4.
Rotate the medium frame so as to setz-axis points down, as shown in Figure 3(a) (red line stands for the axis of the medium frame). Rotate the medium frame with the following rotation angles: 0°, 15°, 30°, 45°, 60°, 90°, 120°, 135°, 150°, 165°, and 180°, and acquire accelerometer outputs from MIMU; denote by [
A
z
x
(
i
)
;
A
z
y
(
i
)
;
A
z
z
(
i
)
].Figure 3
(a)
(b)
(c)Step 5.
Rotate the medium frame and inner frame so as to setx-axis points down, as shown in Figure 3(b). Rotate the medium frame with the rotation angles as in Step 4 and acquire accelerometer outputs from MIMU; denote by [
A
x
x
(
i
)
;
A
x
y
(
i
)
;
A
x
z
(
i
)
].Step 6.
Rotate the medium frame and inner frame so as to sety-axis points down, as shown in Figure 3(c). Rotate the inner frame with the rotation angels as in Step 4 and acquire accelerometer outputs from MIMU; denote by [
A
y
x
(
i
)
;
A
y
y
(
i
)
;
A
y
z
(
i
)
].


### 2.3. Derivation for Parameter Estimation

#### 2.3.1. Derivation for Parameter Estimation in (1)

Parameters in (1) can be estimated by multirate tests. According to (1), when the frame rotates clockwise with the rotation rate i, we can write the following equation: (3)
ω
x
+
i
ω
y
+
i
ω
z
+
i
=
k
g
11
k
g
12
k
g
13
k
g
21
k
g
22
k
g
23
k
g
31
k
g
32
k
g
33
G
x
+
i
G
y
+
i
G
z
+
i
+
ω
x
0
ω
y
0
ω
z
0
.And when the frame rotates counterclockwise with the rotation rate-
i, we can write the following equation:(4)
ω
x
-
i
ω
y
-
i
ω
z
-
i
=
k
g
11
k
g
12
k
g
13
k
g
21
k
g
22
k
g
23
k
g
31
k
g
32
k
g
33
G
x
-
i
G
y
-
i
G
z
-
i
+
ω
x
0
ω
y
0
ω
z
0
.Equations (3)-(4) yield(5)
Δ
ω
x
i
Δ
ω
y
i
Δ
ω
z
i
=
k
g
11
k
g
12
k
g
13
k
g
21
k
g
22
k
g
23
k
g
31
k
g
32
k
g
33
Δ
G
x
i
Δ
G
y
i
Δ
G
z
i
,where Δ
ω
j
(
i
)  (j
=
x
,
y
,
z) means the difference of the rotation rates between the frame rotations clockwise and counterclockwise. Δ
G
i  (
i
=
x
,
y
,
z
) means the difference of gyroscope outputs between the frame rotations clockwise and counterclockwise.According to all acquired data from multirate tests, (5) can be expressed by(6)
Δ
Ω
3
n
×
1
=
Δ
G
3
n
×
9
·
K
G
9
×
1
n
=
29
,where vector K
G
9
×
1 contains the parameters to be calibrated. K
G
9
×
1
=
k
g
11
k
g
12
k
g
13
k
g
21
k
g
22
k
g
23
k
g
31
k
g
32
k
g
33
T. And Δ
Ω
3
n
×
1
=
Δ
Ω
x
Δ
Ω
y
Δ
Ω
z
T. It contains the rotation rates vector when x-axis, y-axis, and z-axis point up. Δ
Ω
x
=
2
×
[
Δ
Ω
x
(
1
)
Δ
Ω
x
(
2
)
Δ
Ω
x
(
3
)
Δ
Ω
x
(
4
)
Δ
Ω
x
(
5
)
Δ
Ω
x
(
6
)
Δ
Ω
x
(
7
)
Δ
Ω
x
(
8
)
Δ
Ω
x
(
9
)
]. Δ
Ω
y
=
2
×
[
Δ
Ω
y
(
1
)
Δ
Ω
y
(
2
)
Δ
Ω
y
(
3
)
Δ
Ω
y
4
Δ
Ω
y
5
Δ
Ω
y
(
6
)
Δ
Ω
y
(
7
)
Δ
Ω
y
(
8
)
Δ
Ω
y
(
9
)
]. Δ
Ω
z
=
2
×
[
Δ
Ω
z
(
1
)
Δ
Ω
z
(
2
)
Δ
Ω
z
(
3
)
Δ
Ω
z
(
4
)
Δ
Ω
z
(
5
)
Δ
Ω
z
(
6
)
Δ
Ω
z
(
7
)
Δ
Ω
z
(
8
)
Δ
Ω
z
(
9
)
Δ
Ω
z
(
10
)
Δ
Ω
z
(
11
)
]. Δ
Ω
i
(
j
) (i
=
x
,
y
,
z, j
=
1
~
11) is shown in appendix. Matrix Δ
G is (7)
Δ
G
=
Δ
G
x
3
×
1
1
0
3
×
1
0
3
×
1
0
3
×
1
Δ
G
x
3
×
1
1
0
3
×
1
0
3
×
1
0
3
×
1
Δ
G
x
3
×
1
1
⋮
Δ
G
y
3
×
1
1
0
3
×
1
0
3
×
1
0
3
×
1
Δ
G
y
3
×
1
1
0
3
×
1
0
3
×
1
0
3
×
1
Δ
G
y
3
×
1
1
⋮
Δ
G
z
3
×
1
1
0
3
×
1
0
3
×
1
0
3
×
1
Δ
G
z
3
×
1
1
0
3
×
1
0
3
×
1
0
3
×
1
Δ
G
z
3
×
1
1
⋮
,
Δ
G
x
3
×
1
i
=
Δ
G
x
x
i
Δ
G
x
y
i
Δ
G
x
z
i
,
Δ
G
y
3
×
1
i
=
Δ
G
y
x
i
Δ
G
y
y
i
Δ
G
y
z
i
,
Δ
G
z
3
×
1
i
=
Δ
G
z
x
i
Δ
G
z
y
i
Δ
G
z
z
i
.By the optimal estimation algorithm such as the least square method, vectorK
G can be obtained by(8)
K
G
9
×
1
=
i
n
v
Δ
G
3
n
×
9
T
·
Δ
G
3
n
×
9
·
Δ
G
3
n
×
9
T
·
Δ
Ω
3
n
×
1
.According to (1), parameters of ω
x
0, ω
y
0, and ω
z
0 can be written as(9)
ω
0
=
ω
-
K
G
9
×
1
·
G
,where ω
0
=
[
ω
x
0
ω
y
0
ω
z
0
]
T, ω
=
[
ω
x
ω
y
ω
z
]
T, and G
=
[
G
x
G
y
G
z
]
T. When the turntable keeps stationary, ω
=
0
3
×
1; then(10)
ω
0
=
-
K
G
9
×
1
·
G
0
.G
(
0
) stands for three gyroscope outputs from MIMU when the turntable keeps stationary. By the data acquired when the turntable keeps stationary, we can calculate ω
0 by(11)
ω
0
=
-
K
G
9
×
1
·
G
x
0
+
G
y
0
+
G
z
0
3
.


#### 2.3.2. Derivation for Parameter Estimation in (2)

Parameters in (2) can be estimated by multiposition tests. Equation (2) can be written as(12)
F
3
×
1
=
A
3
×
12
·
K
A
12
×
1
,where Vector K
A contains the parameters to be calibrated. K
A
12
×
1
=
[
k
a
11
k
a
12
k
a
13
k
a
21
k
a
22
k
a
23
k
a
31
k
a
32
k
a
33
f
x
0
f
y
0
f
z
0
]
T. And F
3
×
1
=
[
f
x
f
y
f
z
]
T. Consider (13)
A
3
×
12
=
A
x
A
y
A
z
0
0
0
0
0
0
1
0
0
0
0
0
A
x
A
y
A
z
0
0
0
0
1
0
0
0
0
0
0
0
A
x
A
y
A
z
0
0
1
.From the acquired data during multiposition tests, we can establish the following estimation mode for the parameters in (2):(14)
F
m
×
1
=
A
m
×
12
·
K
A
12
×
1
m
=
33
,where  F
m
×
1
=
[
F
x
F
y
F
z
]
T. In Figure 3(a), the specific force along three axes during the tests is defined as  F
z
=
[
F
z
1
F
z
2
F
z
3
F
z
4
F
z
5
F
z
6
F
z
7
F
z
(
8
)
F
z
(
9
)
F
z
(
10
)
F
z
(
11
)
]. In Figure 3(b), the specific force along three axes during the tests is defined as F
x
=
[
F
x
(
1
)
F
x
(
2
)
F
x
(
3
)
F
x
(
4
)
F
x
(
5
)
F
x
(
6
)
F
x
(
7
)
F
x
(
8
)
F
x
(
9
)
F
x
(
10
)
F
x
(
11
)
]. In Figure 3(c), the specific force along three axes during the tests is defined as F
y
=
[
F
y
(
1
)
F
y
(
2
)
F
y
(
3
)
F
y
(
4
)
F
y
(
5
)
F
y
(
6
)
F
y
(
7
)
F
y
(
8
)
F
y
(
9
)
F
y
(
10
)
F
y
(
11
)
]. F
i
(
j
) (i
=
x
,
y
,
z, j
=
1
~
11) is shown in appendix. Matrix A
m
×
12 is(15)
A
m
×
12
=
A
x
3
×
1
1
0
3
×
1
0
3
×
1
1
0
0
0
3
×
1
A
x
3
×
1
1
0
3
×
1
0
1
0
0
3
×
1
0
3
×
1
A
x
3
×
1
1
0
0
1
⋮
A
y
3
×
1
1
0
3
×
1
0
3
×
1
1
0
0
0
3
×
1
A
y
3
×
1
1
0
3
×
1
0
1
0
0
3
×
1
0
3
×
1
A
y
3
×
1
1
0
0
1
⋮
A
z
3
×
1
1
0
3
×
1
0
3
×
1
1
0
0
0
3
×
1
A
z
3
×
1
1
0
3
×
1
0
1
0
0
3
×
1
0
3
×
1
A
z
3
×
1
1
0
0
1
⋮
,
A
x
3
×
1
i
=
A
x
x
i
A
x
y
i
A
x
z
i
,
A
y
3
×
1
i
=
A
y
x
i
A
y
y
i
A
y
z
i
,
A
z
3
×
1
i
=
A
z
x
i
A
z
y
i
A
z
z
i
.By the optimal estimation algorithm such as the least square method, vectorK
A
12
×
1 can be obtained by(16)
K
A
12
×
1
=
i
n
v
A
m
×
12
T
·
A
m
×
12
·
A
m
×
12
T
·
F
m
×
1
.


## 2.1. Model Establishment for MIMU

In order to have better application for engineering application, we establish the following calibration model for the outputs of MIMU which contains the triaxial angular rates and specific forces. Different from the existing calibration model, the following model can avoid the calculation of the inverse matrix; thus it can alleviate the computing burden for the system.Output model for triaxial angular rates is as follows:(1)
ω
x
ω
y
ω
z
=
k
g
11
k
g
12
k
g
13
k
g
21
k
g
22
k
g
23
k
g
31
k
g
32
k
g
33
G
x
G
y
G
z
+
ω
x
0
ω
y
0
ω
z
0
.Output model for triaxial specific forces is as follows:(2)
f
x
f
y
f
z
=
k
a
11
k
a
12
k
a
13
k
a
21
k
a
22
k
a
23
k
a
31
k
a
32
k
a
33
A
x
A
y
A
z
+
f
x
0
f
y
0
f
z
0
,where ω
i and f
i
(
i
=
x
,
y
,
z
) mean the output angular rate and specific force of x-axis, y-axis, and z-axis, respectively. G
i and A
i
(
i
=
x
,
y
,
z
) mean the measurement of gyroscope and accelerometer; ω
i
0 and f
i
0
(
i
=
x
,
y
,
z
) mean the bias errors of the angular rate and specific force of x-axis, y-axis, and z-axis. k
i
m
n (i
=
g or a, m
,
n
=
1,2
,
3) is the coefficient related to the measurements of three gyroscopes and accelerometers, which needs to be calibrated.


## 2.2. Calibration Procedure for MIMU

Three-axis turntable applied for the calibration is shown in Figure2. Install MIMU on the plane of the inner frame, while x-axis of MIMU is along with the axis of the medium frame in the turntable, y-axis points to the front and is perpendicular to x-axis, and z-axis is along with the axis of the inner frame in the turntable. Turn the power on, and preheat MIMU for 10 minutes. And then implement the calibration procedure. Calibration procedure includes multirate tests and multiposition tests.Figure 2
Installation of MIMU on the turntable.


### 2.2.1. Multirate Tests

After the initialization of the turntable, multirate tests will be carried out as follows.Step 1.
Adjust the medium frame in order to setx-axis point up, and acquire gyroscope outputs from MIMU when the turntable keeps stationary; denote by G
x
(
0
)
=
[
G
x
x
(
0
)
;
G
x
y
(
0
)
;
G
x
z
(
0
)
]. Rotate the outer frame with the following rotation rates: ±6°/s, ±12°/s, ±24°/s, ±30°/s, ±60°/s, ±120°/s, ±180°/s, ±240°/s, and ±300°/s, and then acquire gyroscope outputs from MIMU. If the rotation rates are positive, denote by [
G
x
x
+
(
i
)
;
G
x
y
+
(
i
)
;
G
x
z
+
(
i
)
]
(
i
=
1
~
9
). If the rotation rates are negative, denote by [
G
x
x
-
(
i
)
;
G
x
y
-
(
i
)
;
G
x
z
-
(
i
)
]
(
i
=
1
~
9
).Step 2.
Adjust the inner frame to sety-axis point up, and acquire the gyroscope outputs from MIMU when the turntable keeps stationary; denote by G
y
(
0
)
=
[
G
y
x
0
;
G
y
y
(
0
)
;
G
y
z
(
0
)
]. By the same way, rotate the outer frame with the rotation rates as in Step 1 and acquire gyroscope outputs from MIMU. If the rotation rates are positive, denote by [
G
y
x
+
(
i
)
;
G
y
y
+
(
i
)
;
G
y
z
+
(
i
)
]
(
i
=
1
~
9
). If the rotation rates are negative, denote by [
G
y
x
-
(
i
)
;
G
y
y
-
(
i
)
;
G
y
z
-
(
i
)
]
(
i
=
1
~
9
).Step 3.
Adjust the medium frame to setz-axis points up, and acquire gyroscope outputs from MIMU when the turntable keeps stationary; denote by G
z
(
0
)
=
[
G
z
x
(
0
)
;
G
z
y
(
0
)
;
G
z
z
(
0
)
]. Rotate the inner frame with the following rotation rates: 0°/s, ±10°/s, ±40°/s, ±100°/s, ±240°/s, ±400°/s, ±720°/s, ±1000°/s, ±1800°/s, ±2160°/s, ±2880°/s, and ±3600°/s, and then acquire gyroscope outputs from MIMU. If the rotation rates are positive, denote by [
G
z
x
+
(
i
)
;
G
z
y
+
(
i
)
;
G
z
z
+
(
i
)
]
(
i
=
1
~
11
). If the rotation rates are negative, denote by [
G
z
x
-
(
i
)
;
G
z
y
-
(
i
)
;
G
z
z
-
(
i
)
]
(
i
=
1
~
11
).


### 2.2.2. Multiposition Tests

After multirate tests have been accomplished, adjust three axes of the turntable as the initialization.Step 4.
Rotate the medium frame so as to setz-axis points down, as shown in Figure 3(a) (red line stands for the axis of the medium frame). Rotate the medium frame with the following rotation angles: 0°, 15°, 30°, 45°, 60°, 90°, 120°, 135°, 150°, 165°, and 180°, and acquire accelerometer outputs from MIMU; denote by [
A
z
x
(
i
)
;
A
z
y
(
i
)
;
A
z
z
(
i
)
].Figure 3
(a)
(b)
(c)Step 5.
Rotate the medium frame and inner frame so as to setx-axis points down, as shown in Figure 3(b). Rotate the medium frame with the rotation angles as in Step 4 and acquire accelerometer outputs from MIMU; denote by [
A
x
x
(
i
)
;
A
x
y
(
i
)
;
A
x
z
(
i
)
].Step 6.
Rotate the medium frame and inner frame so as to sety-axis points down, as shown in Figure 3(c). Rotate the inner frame with the rotation angels as in Step 4 and acquire accelerometer outputs from MIMU; denote by [
A
y
x
(
i
)
;
A
y
y
(
i
)
;
A
y
z
(
i
)
].


## 2.2.1. Multirate Tests

After the initialization of the turntable, multirate tests will be carried out as follows.Step 1.
Adjust the medium frame in order to setx-axis point up, and acquire gyroscope outputs from MIMU when the turntable keeps stationary; denote by G
x
(
0
)
=
[
G
x
x
(
0
)
;
G
x
y
(
0
)
;
G
x
z
(
0
)
]. Rotate the outer frame with the following rotation rates: ±6°/s, ±12°/s, ±24°/s, ±30°/s, ±60°/s, ±120°/s, ±180°/s, ±240°/s, and ±300°/s, and then acquire gyroscope outputs from MIMU. If the rotation rates are positive, denote by [
G
x
x
+
(
i
)
;
G
x
y
+
(
i
)
;
G
x
z
+
(
i
)
]
(
i
=
1
~
9
). If the rotation rates are negative, denote by [
G
x
x
-
(
i
)
;
G
x
y
-
(
i
)
;
G
x
z
-
(
i
)
]
(
i
=
1
~
9
).Step 2.
Adjust the inner frame to sety-axis point up, and acquire the gyroscope outputs from MIMU when the turntable keeps stationary; denote by G
y
(
0
)
=
[
G
y
x
0
;
G
y
y
(
0
)
;
G
y
z
(
0
)
]. By the same way, rotate the outer frame with the rotation rates as in Step 1 and acquire gyroscope outputs from MIMU. If the rotation rates are positive, denote by [
G
y
x
+
(
i
)
;
G
y
y
+
(
i
)
;
G
y
z
+
(
i
)
]
(
i
=
1
~
9
). If the rotation rates are negative, denote by [
G
y
x
-
(
i
)
;
G
y
y
-
(
i
)
;
G
y
z
-
(
i
)
]
(
i
=
1
~
9
).Step 3.
Adjust the medium frame to setz-axis points up, and acquire gyroscope outputs from MIMU when the turntable keeps stationary; denote by G
z
(
0
)
=
[
G
z
x
(
0
)
;
G
z
y
(
0
)
;
G
z
z
(
0
)
]. Rotate the inner frame with the following rotation rates: 0°/s, ±10°/s, ±40°/s, ±100°/s, ±240°/s, ±400°/s, ±720°/s, ±1000°/s, ±1800°/s, ±2160°/s, ±2880°/s, and ±3600°/s, and then acquire gyroscope outputs from MIMU. If the rotation rates are positive, denote by [
G
z
x
+
(
i
)
;
G
z
y
+
(
i
)
;
G
z
z
+
(
i
)
]
(
i
=
1
~
11
). If the rotation rates are negative, denote by [
G
z
x
-
(
i
)
;
G
z
y
-
(
i
)
;
G
z
z
-
(
i
)
]
(
i
=
1
~
11
).


## 2.2.2. Multiposition Tests

After multirate tests have been accomplished, adjust three axes of the turntable as the initialization.Step 4.
Rotate the medium frame so as to setz-axis points down, as shown in Figure 3(a) (red line stands for the axis of the medium frame). Rotate the medium frame with the following rotation angles: 0°, 15°, 30°, 45°, 60°, 90°, 120°, 135°, 150°, 165°, and 180°, and acquire accelerometer outputs from MIMU; denote by [
A
z
x
(
i
)
;
A
z
y
(
i
)
;
A
z
z
(
i
)
].Figure 3
(a)
(b)
(c)Step 5.
Rotate the medium frame and inner frame so as to setx-axis points down, as shown in Figure 3(b). Rotate the medium frame with the rotation angles as in Step 4 and acquire accelerometer outputs from MIMU; denote by [
A
x
x
(
i
)
;
A
x
y
(
i
)
;
A
x
z
(
i
)
].Step 6.
Rotate the medium frame and inner frame so as to sety-axis points down, as shown in Figure 3(c). Rotate the inner frame with the rotation angels as in Step 4 and acquire accelerometer outputs from MIMU; denote by [
A
y
x
(
i
)
;
A
y
y
(
i
)
;
A
y
z
(
i
)
].


## 2.3. Derivation for Parameter Estimation

### 2.3.1. Derivation for Parameter Estimation in (1)

Parameters in (1) can be estimated by multirate tests. According to (1), when the frame rotates clockwise with the rotation rate i, we can write the following equation: (3)
ω
x
+
i
ω
y
+
i
ω
z
+
i
=
k
g
11
k
g
12
k
g
13
k
g
21
k
g
22
k
g
23
k
g
31
k
g
32
k
g
33
G
x
+
i
G
y
+
i
G
z
+
i
+
ω
x
0
ω
y
0
ω
z
0
.And when the frame rotates counterclockwise with the rotation rate-
i, we can write the following equation:(4)
ω
x
-
i
ω
y
-
i
ω
z
-
i
=
k
g
11
k
g
12
k
g
13
k
g
21
k
g
22
k
g
23
k
g
31
k
g
32
k
g
33
G
x
-
i
G
y
-
i
G
z
-
i
+
ω
x
0
ω
y
0
ω
z
0
.Equations (3)-(4) yield(5)
Δ
ω
x
i
Δ
ω
y
i
Δ
ω
z
i
=
k
g
11
k
g
12
k
g
13
k
g
21
k
g
22
k
g
23
k
g
31
k
g
32
k
g
33
Δ
G
x
i
Δ
G
y
i
Δ
G
z
i
,where Δ
ω
j
(
i
)  (j
=
x
,
y
,
z) means the difference of the rotation rates between the frame rotations clockwise and counterclockwise. Δ
G
i  (
i
=
x
,
y
,
z
) means the difference of gyroscope outputs between the frame rotations clockwise and counterclockwise.According to all acquired data from multirate tests, (5) can be expressed by(6)
Δ
Ω
3
n
×
1
=
Δ
G
3
n
×
9
·
K
G
9
×
1
n
=
29
,where vector K
G
9
×
1 contains the parameters to be calibrated. K
G
9
×
1
=
k
g
11
k
g
12
k
g
13
k
g
21
k
g
22
k
g
23
k
g
31
k
g
32
k
g
33
T. And Δ
Ω
3
n
×
1
=
Δ
Ω
x
Δ
Ω
y
Δ
Ω
z
T. It contains the rotation rates vector when x-axis, y-axis, and z-axis point up. Δ
Ω
x
=
2
×
[
Δ
Ω
x
(
1
)
Δ
Ω
x
(
2
)
Δ
Ω
x
(
3
)
Δ
Ω
x
(
4
)
Δ
Ω
x
(
5
)
Δ
Ω
x
(
6
)
Δ
Ω
x
(
7
)
Δ
Ω
x
(
8
)
Δ
Ω
x
(
9
)
]. Δ
Ω
y
=
2
×
[
Δ
Ω
y
(
1
)
Δ
Ω
y
(
2
)
Δ
Ω
y
(
3
)
Δ
Ω
y
4
Δ
Ω
y
5
Δ
Ω
y
(
6
)
Δ
Ω
y
(
7
)
Δ
Ω
y
(
8
)
Δ
Ω
y
(
9
)
]. Δ
Ω
z
=
2
×
[
Δ
Ω
z
(
1
)
Δ
Ω
z
(
2
)
Δ
Ω
z
(
3
)
Δ
Ω
z
(
4
)
Δ
Ω
z
(
5
)
Δ
Ω
z
(
6
)
Δ
Ω
z
(
7
)
Δ
Ω
z
(
8
)
Δ
Ω
z
(
9
)
Δ
Ω
z
(
10
)
Δ
Ω
z
(
11
)
]. Δ
Ω
i
(
j
) (i
=
x
,
y
,
z, j
=
1
~
11) is shown in appendix. Matrix Δ
G is (7)
Δ
G
=
Δ
G
x
3
×
1
1
0
3
×
1
0
3
×
1
0
3
×
1
Δ
G
x
3
×
1
1
0
3
×
1
0
3
×
1
0
3
×
1
Δ
G
x
3
×
1
1
⋮
Δ
G
y
3
×
1
1
0
3
×
1
0
3
×
1
0
3
×
1
Δ
G
y
3
×
1
1
0
3
×
1
0
3
×
1
0
3
×
1
Δ
G
y
3
×
1
1
⋮
Δ
G
z
3
×
1
1
0
3
×
1
0
3
×
1
0
3
×
1
Δ
G
z
3
×
1
1
0
3
×
1
0
3
×
1
0
3
×
1
Δ
G
z
3
×
1
1
⋮
,
Δ
G
x
3
×
1
i
=
Δ
G
x
x
i
Δ
G
x
y
i
Δ
G
x
z
i
,
Δ
G
y
3
×
1
i
=
Δ
G
y
x
i
Δ
G
y
y
i
Δ
G
y
z
i
,
Δ
G
z
3
×
1
i
=
Δ
G
z
x
i
Δ
G
z
y
i
Δ
G
z
z
i
.By the optimal estimation algorithm such as the least square method, vectorK
G can be obtained by(8)
K
G
9
×
1
=
i
n
v
Δ
G
3
n
×
9
T
·
Δ
G
3
n
×
9
·
Δ
G
3
n
×
9
T
·
Δ
Ω
3
n
×
1
.According to (1), parameters of ω
x
0, ω
y
0, and ω
z
0 can be written as(9)
ω
0
=
ω
-
K
G
9
×
1
·
G
,where ω
0
=
[
ω
x
0
ω
y
0
ω
z
0
]
T, ω
=
[
ω
x
ω
y
ω
z
]
T, and G
=
[
G
x
G
y
G
z
]
T. When the turntable keeps stationary, ω
=
0
3
×
1; then(10)
ω
0
=
-
K
G
9
×
1
·
G
0
.G
(
0
) stands for three gyroscope outputs from MIMU when the turntable keeps stationary. By the data acquired when the turntable keeps stationary, we can calculate ω
0 by(11)
ω
0
=
-
K
G
9
×
1
·
G
x
0
+
G
y
0
+
G
z
0
3
.


### 2.3.2. Derivation for Parameter Estimation in (2)

Parameters in (2) can be estimated by multiposition tests. Equation (2) can be written as(12)
F
3
×
1
=
A
3
×
12
·
K
A
12
×
1
,where Vector K
A contains the parameters to be calibrated. K
A
12
×
1
=
[
k
a
11
k
a
12
k
a
13
k
a
21
k
a
22
k
a
23
k
a
31
k
a
32
k
a
33
f
x
0
f
y
0
f
z
0
]
T. And F
3
×
1
=
[
f
x
f
y
f
z
]
T. Consider (13)
A
3
×
12
=
A
x
A
y
A
z
0
0
0
0
0
0
1
0
0
0
0
0
A
x
A
y
A
z
0
0
0
0
1
0
0
0
0
0
0
0
A
x
A
y
A
z
0
0
1
.From the acquired data during multiposition tests, we can establish the following estimation mode for the parameters in (2):(14)
F
m
×
1
=
A
m
×
12
·
K
A
12
×
1
m
=
33
,where  F
m
×
1
=
[
F
x
F
y
F
z
]
T. In Figure 3(a), the specific force along three axes during the tests is defined as  F
z
=
[
F
z
1
F
z
2
F
z
3
F
z
4
F
z
5
F
z
6
F
z
7
F
z
(
8
)
F
z
(
9
)
F
z
(
10
)
F
z
(
11
)
]. In Figure 3(b), the specific force along three axes during the tests is defined as F
x
=
[
F
x
(
1
)
F
x
(
2
)
F
x
(
3
)
F
x
(
4
)
F
x
(
5
)
F
x
(
6
)
F
x
(
7
)
F
x
(
8
)
F
x
(
9
)
F
x
(
10
)
F
x
(
11
)
]. In Figure 3(c), the specific force along three axes during the tests is defined as F
y
=
[
F
y
(
1
)
F
y
(
2
)
F
y
(
3
)
F
y
(
4
)
F
y
(
5
)
F
y
(
6
)
F
y
(
7
)
F
y
(
8
)
F
y
(
9
)
F
y
(
10
)
F
y
(
11
)
]. F
i
(
j
) (i
=
x
,
y
,
z, j
=
1
~
11) is shown in appendix. Matrix A
m
×
12 is(15)
A
m
×
12
=
A
x
3
×
1
1
0
3
×
1
0
3
×
1
1
0
0
0
3
×
1
A
x
3
×
1
1
0
3
×
1
0
1
0
0
3
×
1
0
3
×
1
A
x
3
×
1
1
0
0
1
⋮
A
y
3
×
1
1
0
3
×
1
0
3
×
1
1
0
0
0
3
×
1
A
y
3
×
1
1
0
3
×
1
0
1
0
0
3
×
1
0
3
×
1
A
y
3
×
1
1
0
0
1
⋮
A
z
3
×
1
1
0
3
×
1
0
3
×
1
1
0
0
0
3
×
1
A
z
3
×
1
1
0
3
×
1
0
1
0
0
3
×
1
0
3
×
1
A
z
3
×
1
1
0
0
1
⋮
,
A
x
3
×
1
i
=
A
x
x
i
A
x
y
i
A
x
z
i
,
A
y
3
×
1
i
=
A
y
x
i
A
y
y
i
A
y
z
i
,
A
z
3
×
1
i
=
A
z
x
i
A
z
y
i
A
z
z
i
.By the optimal estimation algorithm such as the least square method, vectorK
A
12
×
1 can be obtained by(16)
K
A
12
×
1
=
i
n
v
A
m
×
12
T
·
A
m
×
12
·
A
m
×
12
T
·
F
m
×
1
.


## 2.3.1. Derivation for Parameter Estimation in (1)

Parameters in (1) can be estimated by multirate tests. According to (1), when the frame rotates clockwise with the rotation rate i, we can write the following equation: (3)
ω
x
+
i
ω
y
+
i
ω
z
+
i
=
k
g
11
k
g
12
k
g
13
k
g
21
k
g
22
k
g
23
k
g
31
k
g
32
k
g
33
G
x
+
i
G
y
+
i
G
z
+
i
+
ω
x
0
ω
y
0
ω
z
0
.And when the frame rotates counterclockwise with the rotation rate-
i, we can write the following equation:(4)
ω
x
-
i
ω
y
-
i
ω
z
-
i
=
k
g
11
k
g
12
k
g
13
k
g
21
k
g
22
k
g
23
k
g
31
k
g
32
k
g
33
G
x
-
i
G
y
-
i
G
z
-
i
+
ω
x
0
ω
y
0
ω
z
0
.Equations (3)-(4) yield(5)
Δ
ω
x
i
Δ
ω
y
i
Δ
ω
z
i
=
k
g
11
k
g
12
k
g
13
k
g
21
k
g
22
k
g
23
k
g
31
k
g
32
k
g
33
Δ
G
x
i
Δ
G
y
i
Δ
G
z
i
,where Δ
ω
j
(
i
)  (j
=
x
,
y
,
z) means the difference of the rotation rates between the frame rotations clockwise and counterclockwise. Δ
G
i  (
i
=
x
,
y
,
z
) means the difference of gyroscope outputs between the frame rotations clockwise and counterclockwise.According to all acquired data from multirate tests, (5) can be expressed by(6)
Δ
Ω
3
n
×
1
=
Δ
G
3
n
×
9
·
K
G
9
×
1
n
=
29
,where vector K
G
9
×
1 contains the parameters to be calibrated. K
G
9
×
1
=
k
g
11
k
g
12
k
g
13
k
g
21
k
g
22
k
g
23
k
g
31
k
g
32
k
g
33
T. And Δ
Ω
3
n
×
1
=
Δ
Ω
x
Δ
Ω
y
Δ
Ω
z
T. It contains the rotation rates vector when x-axis, y-axis, and z-axis point up. Δ
Ω
x
=
2
×
[
Δ
Ω
x
(
1
)
Δ
Ω
x
(
2
)
Δ
Ω
x
(
3
)
Δ
Ω
x
(
4
)
Δ
Ω
x
(
5
)
Δ
Ω
x
(
6
)
Δ
Ω
x
(
7
)
Δ
Ω
x
(
8
)
Δ
Ω
x
(
9
)
]. Δ
Ω
y
=
2
×
[
Δ
Ω
y
(
1
)
Δ
Ω
y
(
2
)
Δ
Ω
y
(
3
)
Δ
Ω
y
4
Δ
Ω
y
5
Δ
Ω
y
(
6
)
Δ
Ω
y
(
7
)
Δ
Ω
y
(
8
)
Δ
Ω
y
(
9
)
]. Δ
Ω
z
=
2
×
[
Δ
Ω
z
(
1
)
Δ
Ω
z
(
2
)
Δ
Ω
z
(
3
)
Δ
Ω
z
(
4
)
Δ
Ω
z
(
5
)
Δ
Ω
z
(
6
)
Δ
Ω
z
(
7
)
Δ
Ω
z
(
8
)
Δ
Ω
z
(
9
)
Δ
Ω
z
(
10
)
Δ
Ω
z
(
11
)
]. Δ
Ω
i
(
j
) (i
=
x
,
y
,
z, j
=
1
~
11) is shown in appendix. Matrix Δ
G is (7)
Δ
G
=
Δ
G
x
3
×
1
1
0
3
×
1
0
3
×
1
0
3
×
1
Δ
G
x
3
×
1
1
0
3
×
1
0
3
×
1
0
3
×
1
Δ
G
x
3
×
1
1
⋮
Δ
G
y
3
×
1
1
0
3
×
1
0
3
×
1
0
3
×
1
Δ
G
y
3
×
1
1
0
3
×
1
0
3
×
1
0
3
×
1
Δ
G
y
3
×
1
1
⋮
Δ
G
z
3
×
1
1
0
3
×
1
0
3
×
1
0
3
×
1
Δ
G
z
3
×
1
1
0
3
×
1
0
3
×
1
0
3
×
1
Δ
G
z
3
×
1
1
⋮
,
Δ
G
x
3
×
1
i
=
Δ
G
x
x
i
Δ
G
x
y
i
Δ
G
x
z
i
,
Δ
G
y
3
×
1
i
=
Δ
G
y
x
i
Δ
G
y
y
i
Δ
G
y
z
i
,
Δ
G
z
3
×
1
i
=
Δ
G
z
x
i
Δ
G
z
y
i
Δ
G
z
z
i
.By the optimal estimation algorithm such as the least square method, vectorK
G can be obtained by(8)
K
G
9
×
1
=
i
n
v
Δ
G
3
n
×
9
T
·
Δ
G
3
n
×
9
·
Δ
G
3
n
×
9
T
·
Δ
Ω
3
n
×
1
.According to (1), parameters of ω
x
0, ω
y
0, and ω
z
0 can be written as(9)
ω
0
=
ω
-
K
G
9
×
1
·
G
,where ω
0
=
[
ω
x
0
ω
y
0
ω
z
0
]
T, ω
=
[
ω
x
ω
y
ω
z
]
T, and G
=
[
G
x
G
y
G
z
]
T. When the turntable keeps stationary, ω
=
0
3
×
1; then(10)
ω
0
=
-
K
G
9
×
1
·
G
0
.G
(
0
) stands for three gyroscope outputs from MIMU when the turntable keeps stationary. By the data acquired when the turntable keeps stationary, we can calculate ω
0 by(11)
ω
0
=
-
K
G
9
×
1
·
G
x
0
+
G
y
0
+
G
z
0
3
.


## 2.3.2. Derivation for Parameter Estimation in (2)

Parameters in (2) can be estimated by multiposition tests. Equation (2) can be written as(12)
F
3
×
1
=
A
3
×
12
·
K
A
12
×
1
,where Vector K
A contains the parameters to be calibrated. K
A
12
×
1
=
[
k
a
11
k
a
12
k
a
13
k
a
21
k
a
22
k
a
23
k
a
31
k
a
32
k
a
33
f
x
0
f
y
0
f
z
0
]
T. And F
3
×
1
=
[
f
x
f
y
f
z
]
T. Consider (13)
A
3
×
12
=
A
x
A
y
A
z
0
0
0
0
0
0
1
0
0
0
0
0
A
x
A
y
A
z
0
0
0
0
1
0
0
0
0
0
0
0
A
x
A
y
A
z
0
0
1
.From the acquired data during multiposition tests, we can establish the following estimation mode for the parameters in (2):(14)
F
m
×
1
=
A
m
×
12
·
K
A
12
×
1
m
=
33
,where  F
m
×
1
=
[
F
x
F
y
F
z
]
T. In Figure 3(a), the specific force along three axes during the tests is defined as  F
z
=
[
F
z
1
F
z
2
F
z
3
F
z
4
F
z
5
F
z
6
F
z
7
F
z
(
8
)
F
z
(
9
)
F
z
(
10
)
F
z
(
11
)
]. In Figure 3(b), the specific force along three axes during the tests is defined as F
x
=
[
F
x
(
1
)
F
x
(
2
)
F
x
(
3
)
F
x
(
4
)
F
x
(
5
)
F
x
(
6
)
F
x
(
7
)
F
x
(
8
)
F
x
(
9
)
F
x
(
10
)
F
x
(
11
)
]. In Figure 3(c), the specific force along three axes during the tests is defined as F
y
=
[
F
y
(
1
)
F
y
(
2
)
F
y
(
3
)
F
y
(
4
)
F
y
(
5
)
F
y
(
6
)
F
y
(
7
)
F
y
(
8
)
F
y
(
9
)
F
y
(
10
)
F
y
(
11
)
]. F
i
(
j
) (i
=
x
,
y
,
z, j
=
1
~
11) is shown in appendix. Matrix A
m
×
12 is(15)
A
m
×
12
=
A
x
3
×
1
1
0
3
×
1
0
3
×
1
1
0
0
0
3
×
1
A
x
3
×
1
1
0
3
×
1
0
1
0
0
3
×
1
0
3
×
1
A
x
3
×
1
1
0
0
1
⋮
A
y
3
×
1
1
0
3
×
1
0
3
×
1
1
0
0
0
3
×
1
A
y
3
×
1
1
0
3
×
1
0
1
0
0
3
×
1
0
3
×
1
A
y
3
×
1
1
0
0
1
⋮
A
z
3
×
1
1
0
3
×
1
0
3
×
1
1
0
0
0
3
×
1
A
z
3
×
1
1
0
3
×
1
0
1
0
0
3
×
1
0
3
×
1
A
z
3
×
1
1
0
0
1
⋮
,
A
x
3
×
1
i
=
A
x
x
i
A
x
y
i
A
x
z
i
,
A
y
3
×
1
i
=
A
y
x
i
A
y
y
i
A
y
z
i
,
A
z
3
×
1
i
=
A
z
x
i
A
z
y
i
A
z
z
i
.By the optimal estimation algorithm such as the least square method, vectorK
A
12
×
1 can be obtained by(16)
K
A
12
×
1
=
i
n
v
A
m
×
12
T
·
A
m
×
12
·
A
m
×
12
T
·
F
m
×
1
.


## 3. Experiment and Error Analysis

### 3.1. Experiment

Implement the experiment as described in Section2. And repeat the calibration procedure for 5 times to ensure the reliability of the calibration.Calculate the average of each acquired data during multirate tests, and then the parameters in the output model for triaxial angular rates can be obtained according to Section2.3.1, which are listed in Table 1.Table 1
Estimation for parameters in the output mode for triaxial angular rates.
Parameters
Number 1
Number 2
Number 3
Number 4
Number 5
Average
k
g
11
0.99974
0.99973
0.99973
0.99970
0.99969
0.99972
k
g
1
2
−0.01196
−0.01191
−0.01191
−0.01208
−0.01209
−0.01199
k
g
1
3
−0.00059
−0.00056
−0.00063
−0.00058
−0.00058
−0.00059
k
g
21
0.02402
0.02390
0.02394
0.02391
0.02393
0.02394
k
g
22
1.00045
1.00047
1.00044
1.00046
1.00050
1.00046
k
g
23
0.01051
0.01053
0.01052
0.01055
0.01055
0.01053
k
g
31
0.01089
0.01087
0.01075
0.01090
0.01093
0.01087
k
g
32
0.01486
0.01487
0.01487
0.01495
0.01496
0.01490
k
g
33
1.00546
1.00551
1.00590
1.00607
1.00600
1.00579
ω
x
0 (rad/s)
−0.03732
−0.04505
−0.04289
−0.04636
−0.04424
−0.04317
ω
y
0 (rad/s)
−0.14450
−0.14354
−0.13033
−0.11033
−0.09103
−0.12395
ω
z
0 (rad/s)
−0.29703
−0.41339
−0.30428
−0.36567
−0.37728
−0.35153Calculate the average of each acquired data during multiposition tests, and then the parameters in the output model for triaxial specific forces can be obtained according to Section2.3.2, which are listed in Table 2.Table 2
Estimation for parameters in the output for triaxial specific forces.
Parameters
Number 1
Number 2
Number 3
Number 4
Number 5
Average
k
a
11
0.96698
0.96695
0.96676
0.96642
0.96625
0.96667
k
a
12
0.01448
0.01419
0.01444
0.01441
0.01417
0.01434
k
a
13
−0.00145
−0.00245
−0.00191
−0.00096
−0.00221
−0.00180
k
a
21
−0.01760
−0.01753
−0.01748
−0.01769
−0.01769
−0.01760
k
a
22
0.96947
0.96966
0.96998
0.96948
0.96964
0.96965
k
a
23
−0.00835
−0.00860
−0.00852
−0.00859
−0.00851
−0.00851
k
a
31
0.01555
0.014672
0.01467
0.014361
0.01454
0.01476
k
a
32
−0.01525
−0.01447
−0.01490
−0.01503
−0.01505
−0.01494
k
a
33
0.99333
0.99286
0.99299
0.99357
0.99357
0.99326
f
x
0 (mg)
−0.16403
−0.16128
−0.17054
−0.16464
−0.16037
−0.16417
f
y
0 (mg)
0.11392
0.10270
0.11108
0.10385
0.10464
0.10724
f
z
0 (mg)
0.28710
0.256601
0.27716
0.25291
0.25164
0.26508Calculate the average of each parameter by five groups of experimental data, and substitute them into (1) and (2), and then the actual triaxial angular rates and specific forces along three axes of MIMU can be obtained.


### 3.2. Error Analysis

Since all parameters have been estimated, the outputs of triaxial angular rates and specific forces along three axes of MIMU can be compensated according to (1) and (2).To evaluate the performance of the proposed calibration method, the verification experiment is implemented. Take the example for the input angular rates on three axes; the output angular rates on each axis will be illustrated. The test is carried out by rotatingx-axis, y-axis, and z-axis of the turntable with different angular rates. And then acquire the original output of MIMU. Compensate the original acquired data by the proposed method and the traditional method, and the results will be illustrated in Figures 4–6.Figure 4
Comparison of the outputs when the input angular rate is alongx-axis.
(a)
Output angular rate ofx-axis
(b)
Output angular rate ofy-axis
(c)
Output angular rate ofz-axisFigure 5
Comparison of the outputs when the input angular rate is alongy-axis.
(a)
Output angular rate ofx-axis
(b)
Output angular rate ofy-axis
(c)
Output angular rate ofz-axisFigure 6
Comparison of the outputs when the input angular rate is alongz-axis.
(a)
Output angular rate ofx-axis
(b)
Output angular rate ofy-axis
(c)
Output angular rate ofz-axisIn Figure4, the input angular rate is along x-axis. The range is −300°/s~300°/s. And the output angular rates on each axis are shown in Figure 4.In Figure5, the input angular rate is along y-axis. The range is −300°/s~300°/s. And the output angular rates on each axis are shown in Figure 5.In Figure6, the input angular rate is along z-axis. The range is −3600°/s~3600°/s. And the output angular rates on each axis are shown in Figure 6.Output errors of the angular rates before and after compensation by the proposed method and traditional method are listed in Table3.Table 3
Output errors of the angular rates.
Input axis
Output error of the angular rate (°/s)
x-axis
y-axis
z-axis
x
Uncompensated
0.004~0.288
−7.305~−7.185
−3.474~2.825
Proposed method
−0.142~0.144
−0.005~0.034
−0.081~0.039
Traditional method
−0.142~0.143
−0.0048~0.034
−0.081~0.041
y
Uncompensated
−3.445~3.738
−0.364~0.4684
−4.829~4.105
Proposed method
−0.005~0.014
−0.15~0.162
0.007~0.069
Traditional method
−0.005~0.015
−0.15~0.161
0.007~0.071
z
Uncompensated
−37.773~37.959
−1.174~1.607
−17.449~16.48
Proposed method
−0.149~0.218
−0.214~0.358
−2.954~4.838
Traditional method
−0.149~0.218
−0.214~0.358
−2.950~4.835In Figure4, the output angular rates of y-axis and z-axis after compensation will not increase as the input angular rate along x-axis increases. The same effect can be illustrated in Figures 5 and 6. In Table 3, it is obvious that the output errors of the angular rate have been well compensated, and the proposed method has the same compensation effect as the traditional method. However, we should emphasize the fact that the proposed method can avoid the calculation of the inverse matrix in the traditional method; thus it can alleviate the computing burden for the system.


## 3.1. Experiment

Implement the experiment as described in Section2. And repeat the calibration procedure for 5 times to ensure the reliability of the calibration.Calculate the average of each acquired data during multirate tests, and then the parameters in the output model for triaxial angular rates can be obtained according to Section2.3.1, which are listed in Table 1.Table 1
Estimation for parameters in the output mode for triaxial angular rates.
Parameters
Number 1
Number 2
Number 3
Number 4
Number 5
Average
k
g
11
0.99974
0.99973
0.99973
0.99970
0.99969
0.99972
k
g
1
2
−0.01196
−0.01191
−0.01191
−0.01208
−0.01209
−0.01199
k
g
1
3
−0.00059
−0.00056
−0.00063
−0.00058
−0.00058
−0.00059
k
g
21
0.02402
0.02390
0.02394
0.02391
0.02393
0.02394
k
g
22
1.00045
1.00047
1.00044
1.00046
1.00050
1.00046
k
g
23
0.01051
0.01053
0.01052
0.01055
0.01055
0.01053
k
g
31
0.01089
0.01087
0.01075
0.01090
0.01093
0.01087
k
g
32
0.01486
0.01487
0.01487
0.01495
0.01496
0.01490
k
g
33
1.00546
1.00551
1.00590
1.00607
1.00600
1.00579
ω
x
0 (rad/s)
−0.03732
−0.04505
−0.04289
−0.04636
−0.04424
−0.04317
ω
y
0 (rad/s)
−0.14450
−0.14354
−0.13033
−0.11033
−0.09103
−0.12395
ω
z
0 (rad/s)
−0.29703
−0.41339
−0.30428
−0.36567
−0.37728
−0.35153Calculate the average of each acquired data during multiposition tests, and then the parameters in the output model for triaxial specific forces can be obtained according to Section2.3.2, which are listed in Table 2.Table 2
Estimation for parameters in the output for triaxial specific forces.
Parameters
Number 1
Number 2
Number 3
Number 4
Number 5
Average
k
a
11
0.96698
0.96695
0.96676
0.96642
0.96625
0.96667
k
a
12
0.01448
0.01419
0.01444
0.01441
0.01417
0.01434
k
a
13
−0.00145
−0.00245
−0.00191
−0.00096
−0.00221
−0.00180
k
a
21
−0.01760
−0.01753
−0.01748
−0.01769
−0.01769
−0.01760
k
a
22
0.96947
0.96966
0.96998
0.96948
0.96964
0.96965
k
a
23
−0.00835
−0.00860
−0.00852
−0.00859
−0.00851
−0.00851
k
a
31
0.01555
0.014672
0.01467
0.014361
0.01454
0.01476
k
a
32
−0.01525
−0.01447
−0.01490
−0.01503
−0.01505
−0.01494
k
a
33
0.99333
0.99286
0.99299
0.99357
0.99357
0.99326
f
x
0 (mg)
−0.16403
−0.16128
−0.17054
−0.16464
−0.16037
−0.16417
f
y
0 (mg)
0.11392
0.10270
0.11108
0.10385
0.10464
0.10724
f
z
0 (mg)
0.28710
0.256601
0.27716
0.25291
0.25164
0.26508Calculate the average of each parameter by five groups of experimental data, and substitute them into (1) and (2), and then the actual triaxial angular rates and specific forces along three axes of MIMU can be obtained.


## 3.2. Error Analysis

Since all parameters have been estimated, the outputs of triaxial angular rates and specific forces along three axes of MIMU can be compensated according to (1) and (2).To evaluate the performance of the proposed calibration method, the verification experiment is implemented. Take the example for the input angular rates on three axes; the output angular rates on each axis will be illustrated. The test is carried out by rotatingx-axis, y-axis, and z-axis of the turntable with different angular rates. And then acquire the original output of MIMU. Compensate the original acquired data by the proposed method and the traditional method, and the results will be illustrated in Figures 4–6.Figure 4
Comparison of the outputs when the input angular rate is alongx-axis.
(a)
Output angular rate ofx-axis
(b)
Output angular rate ofy-axis
(c)
Output angular rate ofz-axisFigure 5
Comparison of the outputs when the input angular rate is alongy-axis.
(a)
Output angular rate ofx-axis
(b)
Output angular rate ofy-axis
(c)
Output angular rate ofz-axisFigure 6
Comparison of the outputs when the input angular rate is alongz-axis.
(a)
Output angular rate ofx-axis
(b)
Output angular rate ofy-axis
(c)
Output angular rate ofz-axisIn Figure4, the input angular rate is along x-axis. The range is −300°/s~300°/s. And the output angular rates on each axis are shown in Figure 4.In Figure5, the input angular rate is along y-axis. The range is −300°/s~300°/s. And the output angular rates on each axis are shown in Figure 5.In Figure6, the input angular rate is along z-axis. The range is −3600°/s~3600°/s. And the output angular rates on each axis are shown in Figure 6.Output errors of the angular rates before and after compensation by the proposed method and traditional method are listed in Table3.Table 3
Output errors of the angular rates.
Input axis
Output error of the angular rate (°/s)
x-axis
y-axis
z-axis
x
Uncompensated
0.004~0.288
−7.305~−7.185
−3.474~2.825
Proposed method
−0.142~0.144
−0.005~0.034
−0.081~0.039
Traditional method
−0.142~0.143
−0.0048~0.034
−0.081~0.041
y
Uncompensated
−3.445~3.738
−0.364~0.4684
−4.829~4.105
Proposed method
−0.005~0.014
−0.15~0.162
0.007~0.069
Traditional method
−0.005~0.015
−0.15~0.161
0.007~0.071
z
Uncompensated
−37.773~37.959
−1.174~1.607
−17.449~16.48
Proposed method
−0.149~0.218
−0.214~0.358
−2.954~4.838
Traditional method
−0.149~0.218
−0.214~0.358
−2.950~4.835In Figure4, the output angular rates of y-axis and z-axis after compensation will not increase as the input angular rate along x-axis increases. The same effect can be illustrated in Figures 5 and 6. In Table 3, it is obvious that the output errors of the angular rate have been well compensated, and the proposed method has the same compensation effect as the traditional method. However, we should emphasize the fact that the proposed method can avoid the calculation of the inverse matrix in the traditional method; thus it can alleviate the computing burden for the system.


## 4. Conclusion

In this paper, we proposed a new calibration method for MIMU which is used for the application of guided spinning projectiles. From the point view of the engineering application, calibration model was established, and the derivation for parameters estimation has been introduced. By multirate tests and multiposition tests, all the parameters in the calibration model were well calibrated. Verification experiment shows that the proposed method has the same compensation effect as the traditional method, but it can alleviate the computing burden for the system. Thus, for the future engineering calibration, the proposed method will have a wide application prospect.



---

*Source: 9108197-2016-06-23.xml*