An Optimal Calibration Method for a MEMS Inertial Measurement Unit

Abstract

An optimal calibration method for a micro-electro-mechanical inertial measurement unit (MIMU) is presented in this paper. The accuracy of the MIMU is highly dependent on calibration to remove the deterministic errors of systematic errors, which also contain random errors. The overlapping Allan variance is applied to characterize the types of random error terms in the measurements. The calibration model includes package misalignment error, sensor-to-sensor misalignment error and bias, and a scale factor is built. The new concept of a calibration method, which includes a calibration scheme and a calibration algorithm, is proposed. The calibration scheme is designed by D-optimal and the calibration algorithm is deduced by a Kalman filter. In addition, the thermal calibration is investigated, as the bias and scale factor varied with temperature. The simulations and real tests verify the effectiveness of the proposed calibration method and show that it is better than the traditional method.

1. Introduction

Nowadays, the progress in micro electro-mechanical systems (MEMSs) has enabled the development of low-cost inertial measurement units (IMUs), which are increasingly used in biomedical applications ^[1], in robots' navigation systems ^[2] and in microsurgical instruments ^[3]. However, the performance of MEMS inertial sensors is degraded by fabrication defects, including produced asymmetric structures, the misalignment of actuation mechanisms and deviations of the centre of mass from the geometric centre. The inner factors cause output errors, which are often divided into deterministic and random errors ^[4]. Moreover, the inertial measurement unit integrated by MEMSs' inertial sensors (MISs) may still contain IMU package misalignment errors and IMU sensor-to-sensor misalignment errors ^[4]. Therefore, it is essential to develop an effective calibration method to reduce the errors and increase the MIMU's precision and stability.

The method most commonly used for the calibration of MEMS accelerometers is the six-position method ^[5], which requires the inertial system to be mounted on a levelled surface with each sensitivity axis of each sensor pointing, alternately, up and down. This calibration method can be used to determine the bias and scale factors of the sensors, but cannot estimate the axes misalignments or non-orthogonalities. The multi-position calibration method is proposed to detect more errors ^[6]. These methods depend on the earth's gravity as a stable physical calibration standard. Furthermore, some special apparatus, such as motion sensing equipment ^[7] and robotic arms ^[8–9], are used for calibration. For the calibration of low-cost MEMS gyroscopes, the Earth's rotation rate is smaller than its resolution. Therefore, the traditional calibration method is dependent on a mechanical platform, rotating the IMU into different, predefined, precisely-controlled orientations and angular rates. Such a method was primarily designed for in-lab tests, which often require the use of expensive equipment ^[5]. Therefore, some calibration methods that do not require the mechanical platform have been proposed. In [10], an optical tracking system is used. In [11], an affordable three-axis rotating platform is designed for the calibration. Meanwhile, schemes for in-field user calibration without external equipment are proposed. Fong et al. ^[12] calibrated gyroscopes by comparing the outputs of the accelerometer and the IMU orientation integration algorithm after arbitrary motions, which requires an initial rough estimate of the gyroscope's parameters. Jurman ^[13] and Hwangbo ^[14] proposed a kind of shape-from-motion calibration method with magnitude constraint of motion. Furthermore, the calculation algorithm is another important issue. The calibration parameters are computed by the algorithm. The least squares method is the algorithm most commonly used in scalar calibration to estimate the calibration parameters ^{[6, 11–14]}. H.L Zhang et al. ^[15] implemented an optimal calibration scheme by maximizing the sensitivity of the measurement norms with respect to the calibration parameters. The algorithms typically lead to a biased estimate of the calibration parameters and may give non-optimal estimates of the calibration coefficients. To avoid this, Panahandeh et al. ^[16] solved the identification problem by using the maximum likelihood estimation (MLE) framework, but only simulation results are presented. In fact, the biases and scale factors vary with temperature. The thermal calibration is also an indispensable process for MIMU. However, the literature seldom completes this type of calibration.

In this paper, an optimal calibration method for a MIMU is presented. Firstly, the measurement noise of the sensors is analysed to provide the information for the calibration. Next, the concept of the calibration method is introduced, which consists of a calibration scheme and a calibration algorithm. The optimal calibration scheme is designed by the D-optimal method ^[17–18]. Meanwhile, the optimal calibration algorithm is deduced by a Kalman filter. The calibration parameters, like scale factors, misalignments and biases, are estimated. Afterwards, thermal calibration is implemented to determine the scale factors and biases in different temperatures. In the end, the results of the simulations and experiments verify the effectiveness of the proposed method, and the performance of the MIMU is improved after the calibration.

The paper is organized as follows. Section 2 presents the hardware of the MIMU and introduces the inertial sensor model. Section 3 describes the character of the sensors' measurement noise, and the overlapping Allan variance is applied. Section 4 gives the MIMU model, where the scale factors, package misalignments, sensor-to-sensor misalignments and biases of the accelerometer and gyroscope triads are considered as the calibration parameters. Then, the optimal calibration algorithm for the gyroscope triad and the accelerometer triad is proposed, followed by the optimal calibration scheme designed by the D-optimal method. The thermal calibration of the MIMU is presented at the end of the section. Section 5 reports the calibration results of MIMU through both simulations and real tests. It shows that the improved multi-position approach outperforms the traditional calibration method. Conclusions are drawn in section 6.

2. Hardware description

A MIMU has been constructed by using MEMS gyros and a MEMS accelerometer. A MMA7260 accelerometer measures the triple-axis acceleration; three single-axis ADXRS300 gyros measure the angular rate in yaw, roll and pitch. Finally, the printed circuit boards (PCBs) of the x, y gyros are kept orthogonal by a slot.

The dynamic range of the ADXRS300 gyro is ±300°/s ^[19], and full range-scale of the MMA7260QT accelerometer is ±1.5 g ^[20]. The outputs from the sensors are analogue voltages that are proportional to the inputs as the rotation rate or acceleration. At zero input, the nominal output of the sensor is bias; positive rotation (clockwise) or acceleration increases the voltage from the nominal null offset value, whereas negative rotation (anticlockwise) or deceleration decreases the voltage output. The null offset value must be subtracted from the raw sensor output measurements. The specifications of the sensors are shown in Table 1.

Table 1.

The specifications of the sensors

	Nominal bias	Nominal Scale Factor
ADXRS300 gyro	2.5 V	5 mv/°/s
MMA7260QT accelerometer	1.65 V	800 mV/g

Based on the above analysis, the output of acceleration measured by the accelerometer can be described by the following equation:

$V_a=K_{a}a+b_a+{\epsilon }_a$

(1)

where V_a is the voltage of acceleration, a is the true acceleration, b_a is the sensor bias, K_a is the scale factor (or acceleration gain) and ε_a is the sensor's noise.

A similar equation can be used to describe the angular velocity measured by the gyros on the single axis:

$V_{\omega }=K_{\omega }\omega +b_{\omega }+{\epsilon }_{\omega }$

(2)

where V _ω is the voltage of angular velocity, ω is the true angular velocity, b_ω is the sensor bias, K_ω is the scale factor (or angular velocity gain) and ε _ω is the sensor's noise.

Although the MEMS gyro's and MEMS accelerometer's parameters are given in the manuals, a dynamic range modification of the MEMS sensors exists. The errors of the parameters as biases and scale factors are deterministic errors, which are reduced by calibration. In addition, the random errors are measurement noises, which should not be ignored in calibration. Therefore, the noise needs to be characterized. For this unit, the data were collected using a multifunction 12-b data-acquisition card from National Instruments, the DAQCard-6008, which is controlled by a script based on National Instruments' LabVIEW 2007.

3. Analysis of the measurement noise

Several variance techniques have been devised for the stochastic modelling of the random errors. These are basically very similar and primarily differ in that various signal processing - by way of weighting functions, window functions, etc. - are incorporated into the analysis algorithms in order to achieve a desired result for improving the model characterizations. The simplest is the Allan variance.

The Allan variance is a method of representing the root mean square (RMS) random drift error as a function of averaging time. It is simple to compute and relatively simple to interpret and understand. The Allan variance method can be used to determine the characteristics of the underlying random processes that give rise to the data noise. This technique can be used to characterize various types of error terms in the inertial-sensor data by performing certain operations on the entire length of data.

Assume that there are N consecutive data points, each having a sample time of t₀. Forming a group of n consecutive data points (with n < N / 2), each member of the group is a cluster. Associated with each cluster is a time T, which is equal to nt₀. If the instantaneous output rate of the inertial sensor is Θ(T), the cluster average is defined as:

${\bar{\mathrm{\Theta }}}_k(T)=\frac{1}{T}{\displaystyle \underset{t_k}{\overset{t_k+T}{\int }}\mathrm{\Theta }(T)dt}$

(3)

where ${\bar{\Theta }}_k(T)$ represents the cluster average of the output rate for a cluster, which starts from the k th data point and contains n data points.

The definition of the subsequent cluster average is:

${\bar{\mathrm{\Theta }}}_{next}\left(T\right)=\frac{1}{T}{\displaystyle \underset{t_{k+1}}{\overset{t_{k+1}+T}{\int }}\mathrm{\Theta }\left(T\right)dt}$

(4)

where t_k+1 = t_k + T.

Performing the average operation for each of the two adjacent clusters can form the difference:

${\mathit{\xi }}_{k+1,k}={\bar{\mathrm{\Theta }}}_{next}\left(T\right)-{\bar{\mathrm{\Theta }}}_k\left(T\right)$

(5)

For each cluster time T, the ensemble of ξ _s defined by (5) forms a set of random variables. The quantity of interest is the variance of ξ _s over all the clusters of the same size that can be formed from the entire data. Thus, the Allan variance of length T is defined as ^[21]:

${\mathit{\sigma }}^2(T)=\frac{1}{2\left(N-2n\right)}{\displaystyle \sum _{k=1}^{N-2n}{\left[{\bar{\mathrm{\Theta }}}_{next}(T)-{\bar{\mathrm{\Theta }}}_k(T)\right]}^2}$

(6)

Obviously, for any finite number of data points (N), a finite number of clusters of a fixed length (T) can be formed. Hence, equation (6) represents an estimation of the quantity σ ² (T) whose quality of estimate depends on the number of independent clusters of a fixed length that can be formed. It is normally plotted as the square root of the Allan variance σ (T) versus T on a log-log plot.

In addition, the percentage error is equal to:

$\mathit{\sigma }(\mathit{\delta })=\frac{1}{\sqrt{2\left(N/n-1\right)}}$

(7)

where N is the total number of data points in the entire run and n is the number of data points contained in the cluster.

Equation (7) shows that the estimation errors in the region of a short cluster length T are small where the number of independent clusters in these regions is large. In contrast, the estimation errors in the region of a long cluster length T are large where the number of independent clusters in these regions is small. In order to avoid the estimation error increasing due to a lack of samples, the overlapping Allan variance was proposed ^[22] The use of overlapping samples improves the confidence of the resulting stability estimate. The formula is expressed as follows:

${\mathit{\sigma }}^2(T)=\frac{1}{2m^2(N-2m+1)}{\displaystyle \sum _{j=1}^{N-2m+1}{\left\{{\displaystyle \sum _{i=j}^{j+m-1}\left({\bar{\mathrm{\Theta }}}_{i+m}\left(T\right)-{\bar{\mathrm{\Theta }}}_i\left(T\right)\right)}\right\}}^2}$

(8)

where m is the averaging factor.

The static data from the MIMU were collected with a sampling rate of 1000 Hz at room temperature 25°C. By applying the overlapping Allan-variance method to the whole data set, a log–log plot of the overlapping Allan standard deviation versus the cluster time is shown in Figure 1 for the accelerometers data and Figure 2 for the gyros data. As shown in these figures, the drift characteristic of the sensors is at its worst as σ takes the longest time interval - of about 20 s - to converge, which means that the bias drifts of the accelerometers and gyros are much noisier. This implies that the sensor bias should be averaged over a period of at least 20 s so that the average bias will not change significantly over the next few 20 s intervals. Within 20 s, the random error can be considered as white noise.

Figure 1.

Accelerometers overlapping Allan variance results

Figure 2.

Gyros overlapping Allan variance results

In dealing with the sensor signals over a period of time determined using the Allan variance analysis above, it will keep the bias drift minimal during the successive time period when the calibration data are collected. Accordingly, the white noise is only considered in the calibration.

4. Calibration method

Calibration is the process of comparing the instruments output with known reference information and determining the coefficients that force the output to agree with the reference information over a range of output values ^[23]. Moreover, many researchers ^[5–16] have developed calibration methods for the IMU. Before the description of our own calibration method, a clear concept of the calibration method is introduced first of all. The calibration method of the IMU should consist of two aspects, which are a calibration scheme and a calibration algorithm. The scheme is to design the experiments, while the algorithm is to compute the parameters by the experimental data. A novel calibration method is presented in this section according to this concept.

Firstly, the complete calibration model of the MIMU is described. Next, an optimal calibration method is proposed by an optimal calibration scheme and an optimal calibration algorithm. The scheme is the design of a calibration experiment for the gyroscope triad and the accelerometer triad, and the calibration results are computed by the optimal algorithm. Finally, thermal calibration is described.

4.1 Calibration model

The measurement model of the accelerometer or gyroscope generally includes bias, scale factors and random error. Bias and scale factors are considered as the most significant parameters, which change with different temperatures. It is also a very important characteristic for inertial sensors (though this is dealt with by independent thermal calibration, as explained at the end). In this section, the calibration procedure is performed at a stable temperature and the self-heating effects can be ignored, as the MIMU is calibrated after the sensors are warmed up to thermal stability. Moreover, the random error can be considered as white Gaussian noise from the last section.

In this paper, the MIMU consists of three almost orthogonally mounted gyroscopes and a three-axis accelerometer. As to the unit, the misalignment error would be induced because of the sensors' installation. The misalignment error is divided into two sources: package misalignment error and sensor-to-sensor misalignment error. Package misalignment error is defined as the angle between the true axis of sensitivity and the body axis of the package. Sensor-to-sensor misalignment error defines the misalignment error due to the non-orthogonality of the IMU's axes.

The package misalignment angle can be defined as three angles. First, for a rotation about the z -axis by an angle of θ _z , the matrix can be represented as follows:

$C_1=\left[\begin{array}{ccc}\cos {\mathit{\theta }}_z& \sin {\mathit{\theta }}_z& 0\\ -\sin {\mathit{\theta }}_z& \cos {\mathit{\theta }}_z& 0\\ 0& 0& 1\end{array}\right]$

(9)

Second, for a rotation about the y -axis by an angle θ _y, the matrix is represented as follows:

$C_2=\left[\begin{array}{ccc}\cos {\mathit{\theta }}_y& 0& -\sin {\mathit{\theta }}_y\\ 0& 1& 0\\ \sin {\mathit{\theta }}_y& 0& \cos {\mathit{\theta }}_y\end{array}\right]$

(10)

And finally, for a rotation about the x -axis by an angle of θ _x , the matrix is represented as follows:

$C_3=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos {\mathit{\theta }}_x& -\sin {\mathit{\theta }}_x\\ 0& \sin {\mathit{\theta }}_x& \cos {\mathit{\theta }}_x\end{array}\right]$

(11)

The three angles can be considered as small angles. Then, the package misalignment error can be represented by the following equation:

$M_P=C_1{C}_2{C}_3\approx \left[\begin{array}{ccc}1& -{\mathit{\theta }}_z& {\mathit{\theta }}_y\\ {\mathit{\theta }}_z& 1& -{\mathit{\theta }}_x\\ -{\mathit{\theta }}_y& {\mathit{\theta }}_x& 1\end{array}\right]$

(12)

The sensor-to-sensor misalignment error or non-orthogonality can be also considered as creating a series of rotation matrices that define the relationship of the misaligned axes to those of the perfectly orthogonal triad. Generally, it is considered that the axes x of the frames coincide and - assuming small angles - the non-orthogonality matrix can be approximated by:

$M_o=\left[\begin{array}{ccc}1& 0& 0\\ {\mathit{\beta }}_{yz}& 1& 0\\ -{\mathit{\beta }}_{zy}& {\mathit{\beta }}_{zx}& 1\end{array}\right]$

(13)

Then, the total misalignment error can be written as:

$\begin{array}{l}M=M_o{M}_p\\ =\left[\begin{array}{ccc}1& -{\mathit{\theta }}_z& {\mathit{\theta }}_y\\ {\mathit{\beta }}_{yz}+{\mathit{\theta }}_z& 1-{\mathit{\beta }}_{yz}{\mathit{\theta }}_z& {\mathit{\beta }}_{yz}{\mathit{\theta }}_y-{\mathit{\theta }}_x\\ -{\mathit{\beta }}_{zy}+{\mathit{\beta }}_{zx}{\mathit{\theta }}_z-{\mathit{\theta }}_y& {\mathit{\beta }}_{yz}{\mathit{\theta }}_z+{\mathit{\beta }}_{zx}+{\mathit{\theta }}_x& -{\mathit{\beta }}_{zy}{\mathit{\theta }}_y-{\mathit{\beta }}_{zx}{\mathit{\theta }}_x+1\end{array}\right]\\ \approx \left[\begin{array}{ccc}1& -{\mathit{\theta }}_z& {\mathit{\theta }}_y\\ {\mathit{\beta }}_{yz}+{\mathit{\theta }}_z& 1& -{\mathit{\theta }}_x\\ -{\mathit{\beta }}_{zy}-{\mathit{\theta }}_y& {\mathit{\beta }}_{zx}+{\mathit{\theta }}_x& 1\end{array}\right]\end{array}$

(14)

Consequently, the measurement of the accelerometer cluster can be expressed as:

$a_c=K_a{M}_a{a}_g+b_a+v_a$

(15)

where a _g denotes the input specific force expressed in platform coordinates and v _a is the measurement noise. The scale factor matrix and bias vector of the accelerometer are defined, respectively, as:

$K_a=diag(k_{xa},k_{ya},k_{za}),b_a={\left[\begin{array}{ccc}{b}_{xa}& b_{ya}& b_{za}\end{array}\right]}^{\mathrm{T}}$

where k_ia and b_ia, respectively, denote the scale factor and the bias of the i th accelerometer output, i = x,y,z. The misalignment matrix of the accelerometer is written as:

${\mathit{M}}_a=\left[\begin{array}{ccc}1& -{\hspace{0.17em}}^a{\mathit{\theta }}_z& {\hspace{0.17em}}^a{\mathit{\theta }}_y\\ {\hspace{0.17em}}^a{\mathit{\beta }}_{yz}+{\hspace{0.17em}}^a{\mathit{\theta }}_z& 1& -{\hspace{0.17em}}^a{\mathit{\theta }}_x\\ -{\hspace{0.17em}}^a{\mathit{\beta }}_{zy}-{\hspace{0.17em}}^a{\mathit{\theta }}_y& {\hspace{0.17em}}^a{\mathit{\beta }}_{zx}+{\hspace{0.17em}}^a{\mathit{\theta }}_x& 1\end{array}\right]$

(16)

Analogously, the measurement of the gyroscope cluster can be written as:

${\mathit{\omega }}_c=K_g{M}_g{\mathit{\omega }}_g+b_g+v_g$

(17)

where ω _c denotes the true platform angular velocity with respect to the inertial coordinates expressed in platform coordinates, K _g is the diagonal scale factor matrix, b _g is the bias vector of the gyroscope cluster, and v _g is the measurement noise.

The misalignment matrix of the gyroscope triad is written as:

${\mathit{M}}_g=\left[\begin{array}{ccc}1& -{\hspace{0.17em}}^g{\mathit{\theta }}_z& {\hspace{0.17em}}^g{\mathit{\theta }}_y\\ {\hspace{0.17em}}^g{\mathit{\beta }}_{yz}+{\hspace{0.17em}}^g{\mathit{\theta }}_z& 1& -{\hspace{0.17em}}^g{\mathit{\theta }}_x\\ -{\hspace{0.17em}}^g{\mathit{\beta }}_{zy}-{\hspace{0.17em}}^g{\mathit{\theta }}_y& {\hspace{0.17em}}^g{\mathit{\beta }}_{zx}+{\hspace{0.17em}}^g{\mathit{\theta }}_x& 1\end{array}\right]$

(18)

The task of MIMU calibration is to estimate such parameters as scale factors, misalignments, biases of the accelerometer and the gyroscopes. Meanwhile, the experiment is implemented at 25°C in the lab. The key technology of the calibration can be divided according to two aspects: the calibration scheme and the calibration algorithm. We will discuss these separately in the following sections.

4.2 Calibration scheme

The calibration scheme is the experiment design for the MIMU. Generally, a multi-position calibration scheme is used for the IMU's calibration. The general principle of the method is to design enough positions to estimate the calibration parameters. At least 12 different equations for determining (that is, at least 12 positions for the MIMU) are required in calibration, as there are 12 unknown parameters in both the accelerometer triad and the gyroscope triad. In order to avoid computational singularity in estimation, more positions are desired to get numerically reliable results in reality. Afterwards, the question as to how the positions are designed in calibration arises. The current literature has seldom discussed in detail how the positions are optimized - in other words, the scheme not only makes all the calibration parameters identifiable, but also maximizes their numerical accuracy. In the current context, the optimal calibration positions are proposed by the D-optimal method.

According to equations (15)–(17), the following equation can be deduced:

$O=PI+b+v$

(19)

where P = K _a M _a or K _g M _g , I = a _g or ω _g denotes the input vector of MIMU, O = a _c or ω _c denotes the vector of the MIMU's measurement, v = v _a or v _g , b = b _a or b _g .

By rearranging equation (19), we can attain the following equations:

$\{\begin{array}{c}{\hspace{0.17em}}^{x}O={\mathit{P}}_1\cdot \mathit{I}+b_x+{\hspace{0.17em}}^{x}v\\ {\hspace{0.17em}}^{y}O={\mathit{P}}_2\cdot \mathit{I}+b_y+{\hspace{0.17em}}^{y}v\\ {\hspace{0.17em}}^{z}O={\mathit{P}}_3\cdot \mathit{I}+b_z+{\hspace{0.17em}}^{z}v\end{array}$

(20)

where $\begin{array}{l}{P}_1={\left[\begin{array}{ccc}{P}_{11}& P_{12}& P_{13}\end{array}\right]}^{\mathrm{T}},\\ P_2={\left[\begin{array}{ccc}{P}_{21}& P_{22}& P_{23}\end{array}\right]}^{\mathrm{T}},\\ P_3={\left[\begin{array}{ccc}{P}_{31}& P_{32}& P_{33}\end{array}\right]}^{\mathrm{T}},\end{array}$ ,

ⁱO denotes the MIMU's measurement of the i axis, i=x,y,z.

Assume Ξ is an experimental procedure, which is used to calibrate the parameters in equation (20). The experimental procedure is composed of n tests, defined as U_k (k ∈ [1,n). Each test U_k corresponds to a test position of Ξ, which is generated by the input vector Γ_k = [ I _k 1], where I _k is the vector at test point k. As a result, every experimental procedure can be mathematically represented as:

$\text{Ξ}=\left\{U_k\left({\mathit{\Gamma }}_k\right)|k\in \left[1,n\right]\right\}$

(21)

with associated outputs ⁱ O.

It is shown that a D-optimal design can be achieved with p ≤ n ≤ p(p + 1)/2, where p is the number of parameters to be estimated. From equation (21), one can see that there are four unknown parameters for each output axis, so the optimal number of measurement positions must exist in [5,13].

Then, the optimal function can be described as follows:

$\begin{array}{l}\text{Find}:\hspace{1em}\text{Ξ}\ast =\left[\begin{array}{cccc}{\mathit{\Gamma }}_1& {\mathit{\Gamma }}_2& \cdots & {\mathit{\Gamma }}_k\end{array}\right]\in {\rm E}\\ \text{With}:\hspace{1em}\hspace{1em}\hspace{1em}\text{Ξ}=\arg \underset{\text{Ξ}\in \mathrm{E}}{\max }\{D_n\left(\text{Ξ}\right)=\det \left(F_n^{\hspace{0.17em}\mathrm{T}}{F}_n/n\right)\}\end{array}$

(22)

where E is the global region of positions and F _n = [Γ₁ Γ₂ ⋯ Γ_k] for the design information matrix of n test positions.

In addition, the procedure of D-optimal design is described as follows:

• 1)

Select n test positions from global candidate positions and calculate the object function D _ni

2)

Select a maximum D _ni and let D _n = max D _ni

3)

n ≤ 12,return to step 1)

4)

Optimal experiment positions Ξ*

Four optimal positions for each ⁱO are acquired; therefore, totally, 12 optimal positions can be attained for calibration. Then, there are the duplicate positions. After removing the duplicate positions, the nine optimal positions (which are listed in Table 2) are used for calibration. It is noticed that the MEMS gyros in the paper are low cost and cannot measure the Earth's rotation rate. Therefore, a rotation rate table is used in the calibration to produce a reference signal. Meanwhile, gravity is also considered as a reference signal.

Table 2.

Optimal positions for calibration

No.	MIMU Postures			Illustration
	x	y	z
1	1	0	0
2	−1	0	0
3	$\sqrt{2}/2$	$\sqrt{2}/2$	0
4	$\sqrt{2}/2$	0	$\sqrt{2}/2$
5	0	1	0
6	0	−1	0
7	0	$\sqrt{2}/2$	$\sqrt{2}/2$
8	0	0	1
9	0	0	−1

4.3 Calibration algorithm

The calibration algorithm is the computation of the calibration parameters from the data, which is collected by the calibration scheme. Because the collected data from the MIMU contain random noise, the optimal estimate algorithm is generally adopted. The least squares algorithm is the most commonly used in estimating the calibration parameters. The algorithms typically lead to a biased estimate of the calibration parameters and may give non-optimal estimates of the calibration coefficients. Consequently, a Kalman filter is designed for the calibration algorithm in what follows.

4.3.1 Kalman filter design

A Kalman filter (KF) eliminates random noise and errors using knowledge about the state-space representation of the system and uncertainties in the process: the measurement noise and the process noise. As such, it is very useful as a signal filter in the case of signals with random disturbances or else when signal errors can be treated as an additional state-space variable in the system. For linear Gaussian systems, the Kalman filter is the optimal minimum mean square error estimator.

The model can be expressed as follows:

$\begin{array}{c}{X}_{i+1}={\Phi X}_i+w_i\\ Z_{i+1}=H\cdot X_{i+1}+{\mathit{\epsilon }}_{i+1}\end{array}$

(23)

where X _i is the state vector at time i and:

Z _i + 1 - the system output (measured signal) at time i,

w _i - the process noise at time i,

ε _i -the measurement noise at time i,

Φ , H - matrices of the state-space representation: Φ - the state; H - the output.

Moreover, the Kalman filter can be derived by the following process:

Time updating:

$\{\begin{array}{l}{P}_{i|i-1}={\Phi }_{i-1}{P}_{i-1}{\Phi }_{i-1}^{\hspace{1em}T}\\ {\hat{X}}_{i|i-1}={\hat{X}}_{i-1}\end{array}$

(24)

Measurement updating:

$\{\begin{array}{l}{K}_i=P_{i|i-1}{H}_i^T{\left(H_i{P}_{i|i-1}{H}_i^T+R_i\right)}^{-1}\\ {\widehat{X}}_i={\widehat{X}}_{i|i-1}+K_i\left(Z_i-H_i{\widehat{X}}_{i|i-1}\right)\\ P_i=\left(I-K_i{H}_i\right)P_{i|i-1}\end{array}$

(25)

According to the optimal positions, a bank of Kalman filters is deduced:

$\{\begin{array}{c}{x}_{i+1}={\hspace{0.17em}}^x\Phi x_i\\ {\hspace{0.17em}}^x{Z}_{i+1}={\hspace{0.17em}}^{x}H{x}_{i+1}+{\hspace{0.17em}}^x{\varepsilon }_{i+1}\end{array}$

(26)

$\{\begin{array}{l}{y}_{i+1}={\hspace{0.17em}}^y\Phi y_i\\ {\hspace{0.17em}}^y{Z}_{i+1}={\hspace{0.17em}}^{y}H{y}_{i+1}+{\hspace{0.17em}}^y{\varepsilon }_{i+1}\end{array}$

(27)

$\{\begin{array}{l}{z}_{i+1}={\hspace{0.17em}}^z\Phi z_i\\ {\hspace{0.17em}}^z{Z}_{i+1}={\hspace{0.17em}}^{z}H{z}_{i+1}+{\hspace{0.17em}}^z{\varepsilon }_{i+1}\end{array}$

(28)

where ^x Φ = I _4×4, ^x X = [b_x P₁₁ P_l2 P₁₃]^T,

$\begin{array}{c}{\hspace{0.17em}}^{x}H=\left[\begin{array}{cccc}1& r& 0& 0\\ 1& -r& 0& 0\\ 1& r/\sqrt{2}& r/\sqrt{2}& 0\\ 1& r/\sqrt{2}& 0& r/\sqrt{2}\end{array}\right],r=g\hspace{0.17em}\text{or}\hspace{0.17em}\omega \hspace{0.17em},\\ {\hspace{0.17em}}^{\mathrm{y}}\Phi =I_{4\times 4},{\hspace{0.17em}}^{y}X={\left[\begin{array}{cccc}{b}_y& P_{21}& P_{22}& P_{23}\end{array}\right]}^{\mathrm{T}},\\ {\hspace{0.17em}}^{y}H=\left[\begin{array}{cccc}1& 0& -r& 0\\ 1& 0& -r& 0\\ 1& r/\sqrt{2}& r/\sqrt{2}& 0\\ 1& 0& r/\sqrt{2}& r/\sqrt{2}\end{array}\right],\\ {\hspace{0.17em}}^{\mathrm{z}}\Phi =I_{4\times 4},{\hspace{0.17em}}^{z}X={\left[\begin{array}{cccc}{b}_y& P_{21}& P_{22}& P_{23}\end{array}\right]}^{\mathrm{T}},\\ {\hspace{0.17em}}^{z}H=\left[\begin{array}{cccc}1& 0& 0& r\\ 1& 0& 0& -r\\ 1& 0& r/\sqrt{2}& r/\sqrt{2}\\ 1& r/\sqrt{2}& 0& r/\sqrt{2}\end{array}\right].\end{array}$

Three simpler Kalman filters have been designed for the MIMU calibration. Next, the observability is analysed to make sure that the calibration parameters are identifiable.

According to [24], the n -state discrete linear time-invariant system has the observability matrix G defined by:

$G\left(\mathit{\varphi },H\right)=\left[\begin{array}{cccc}H& H\mathit{\varphi }& \dots & H{\mathit{\varphi }}^{n-1}\end{array}\right]$

(29)

The system is observable if, and only if, rank(G) = n. We can derive the following equations from equations (26)–(28)

$rank\left({\hspace{0.17em}}^{x}H\right)=4,rank\left({\hspace{0.17em}}^{y}H\right)=4,rank\left({\hspace{0.17em}}^{z}H\right)=4$

(30)

Because ⁱ Φ (i = x,y,z) are identity matrices, the rank of the observability matrix is equal to the quantity of states. We can conclude that three Kalman filters are observable - in another words, the parameters can be estimated by the measurements, which are measured by the nine-position experiment.

Therefore, the scale factors, biases and orthogonalization angles of the MIMU are estimated by the three Kalman filters. However, the scale factors and biases are not stable in different temperatures. Therefore, thermal calibration is also important for the MIMU.

4.4 Thermal calibration

The purpose of thermal calibration is to reduce the errors caused by the variation of the scale factors and biases of the MIMU when operated under different temperatures. There are two main approaches for thermal testing ^[5]: (1) Allow the IMU enclosed in the thermal chamber to stabilize at a particular temperature corresponding to the temperature of the thermal chamber and then record the data (this method of recording the data at specific temperature points is called the ‘Soak method’). (2) In the so-called ‘thermal ramp’ method, the IMU temperature is linearly increased or decreased for a certain period of time. We would use the Soak method to investigate the thermal effect of the sensors.

According to equation (15) and equation (17), the model of thermal calibration can be described:

${\mathit{\omega }}_c=K_g^t{M}_g{\mathit{\omega }}_g+b_g^t+v_g$

(31)

$a_c=K_a^t{M}_a{a}_g+b_a^t+v_a$

(32)

where K ^t _g and K ^t _a denote the scale factors at a certain temperature,

and b ^t _g and b ^t _a denote the biases at a certain temperature.

The misalignment parameters can be obtained by the last section and they are not varied by the temperature, so only the scale factors and biases should be estimated. The calibration scheme can be designed by no. 1–6 in Table 2.

Analogously, the calibration algorithms are deduced:

$\{\begin{array}{l}{\hspace{0.17em}}^t{X}_{i+1}^a={\Phi }_a{\hspace{0.17em}}^t{X}_i^a\hfill \\ {\hspace{0.17em}}^t{Z}_{i+1}^a={\hspace{0.17em}}^t{H}_a{\hspace{0.17em}}^t{X}_{i+1}^a+{\hspace{0.17em}}^t{\varepsilon }_{i+1}^a\hfill \end{array}$

(33)

$\{\begin{array}{l}{\hspace{0.17em}}^t{X}_{i+1}^{\omega }={\Phi }_{\omega }{\hspace{0.17em}}^t{X}_i^{\omega }\hfill \\ {\hspace{0.17em}}^t{Z}_{i+1}^{\omega }={\hspace{0.17em}}^t{H}_{\omega }{\hspace{0.17em}}^t{X}_{i+1}^{\omega }+{\hspace{0.17em}}^t{\varepsilon }_{i+1}^{\omega }\hfill \end{array}$

(34)

where ${\hspace{0.17em}}^t{X}^a={[k_x^a{k}_y^a{k}_z^a{b}_x^a{b}_y^a{b}_z^a]}^T$

$\begin{array}{c}{\hspace{0.17em}}^t{X}^{\omega }={\left[\begin{array}{cccccc}{k}_x^{\omega }& k_y^{\omega }& k_z^{\omega }& b_x^{\omega }& b_y^{\omega }& b_z^{\omega }\end{array}\right]}^{\mathrm{T}},\\ {\Phi }^a={\Phi }^{\omega }=I_{6\times 6},\\ {\hspace{0.17em}}^t{H}^a={\hspace{0.17em}}^t{H}^{\omega }=\left[\begin{array}{cccccc}r& 0& 0& 1& 0& 0\\ -r& 0& 0& 1& 0& 0\\ 0& r& 0& 0& 1& 0\\ 0& -r& 0& 0& 1& 0\\ 0& 0& r& 0& 0& 1\\ 0& 0& -r& 0& 0& 1\end{array}\right].\end{array}$

It is noticed that the variances of the MIMU measurements are different for different temperatures, and so the covariance matrices of the measurement noise should be calculated for different temperatures. In thermal calibration, a turntable and a thermal chamber are assembled together to form a thermal-turntable unit, as shown in Figure 3.

Figure 3.

The thermal test setup

5. Calibration Results

5.1 Simulation and calibration results

In order to verify the calibration method, a Monte Carlo simulation is first implemented. The true values of the 12 parameters in the calibration model are listed in Table 3. It is noted that the parameters are dimensionless parameters.

Table 3.

The simulation parameters

	Sensor 1	Sensor 2	Sensor 3
Bias	2.0	2.5	2.9
Scale factor	0.01	0.05	0.09
θ x		0.01
θ y		0.05
θ z		0.1
β yz		0.5
β zy		0.5
β zx		0.8

These parameters used in the simulation were chosen to be within the range of the expected MIMU parameters. The sensor noise standard deviation is:

$\mathit{\sigma }=\left[\begin{array}{cc}0.001& \begin{array}{cc}0.001& 0.001\end{array}\end{array}\right]$

(35)

Both the traditional method and the proposed approach were performed in the simulation. It is noted that the traditional method is the method of least squares ^[5]. The results are shown in Figures 4–7.

Figure 4.

The bias errors of the two calibration methods

Figure 5.

The scale factors error of the two calibration methods

Figure 6.

The package misalignment angle errors of the two calibration methods

Figure 7.

The orthogonalization angle errors of the two calibration methods

The parameters' errors are smaller under the proposed method. Meanwhile, we find that parameters estimated by the traditional method will deviate significantly when the noise standard deviation is increased. However, the proposed method is still effective.

The simulation results have shown the performance of the filter and the accuracy of the proposed method. Next, the MIMU is calibrated in the lab, using a cube and a rate table. The outputs were sampled at a rate of 1000 Hz and the LabVIEW program stored the digitalized analogue measurements from the MIMU. The collected data were later read into the MATLAB program for processing. The calibration results are shown as follows. Figures 8–11 show the parameters of the three-axis accelerometers while the gyros' are shown in Figures 12–15.

Figure 8.

Bias of the accelerometers

Figure 9.

Scale factor of the accelerometers

Figure 10.

Package misalignment angles of the accelerometers

Figure 11.

Non-orthogonality angles of the accelerometers

Figure 12.

Bias of the gyros

Figure 13.

Scale factor of the gyros

Figure 14.

Package misalignment angles of the gyros

Figure 15.

Non-orthogonality angles of the gyros

The estimated results are listed in Table 4. It indicates that the estimating procedure is convergent and it takes less than two seconds to complete. The parameters converge quickly to the stable values.

Table 4.

The calibration results of the MIMU

	Acc. X	Acc. Y	Acc. Z	Gyro X	Gyro Y	Gyro Z
Bias	(m/s2)			(° /s)
	1.79	1.66	1.74	2.51	2.52	2.48
Scale factor	(mv/m/s2)			(mv/° /s)
	80.80	80.14	80.84	5.07	5.09	5.17
θ x (°)	0.06			−0.18
θ y (°)	0.20			0.05
θ z (°)	0.19			0.15
β yz (°)	0.01			1.32
β zy (°)	−0.02			1.05
β zx (°)	−0.01			1.15

5.2 Thermal calibration results

The temperature of the thermal chamber varies from 0°C to 50°C. A total of nine different temperatures are considered in this experiment. For each temperature, the MIMU is allowed to stabilize before recording data.

We calculate the variances of the MIMU measurements at different temperatures. The results are listed in Table 5, and they are different for the different temperatures. The scale factors and biases for the MIMU at different temperatures are shown in Figures 16–19. We can see that the biases and scale factors of the MIMU vary significantly with temperature. Hence, it is necessary to get the thermal calibration model for low-cost MEMS sensors to compensate for these biases and scale factor drift with temperature.

Table 5.

The measurement variance for different temperatures

Temperature(°C)	Variance http://proceedings%.aip%.org/resource/2/apcpcs/(10-5v2)
	Gyro	Gyro	Gyro	Acc.	Acc.	Acc.
	X	Y	Z	X	Y	Z
0	2.14	2.46	1.77	5.92	4.86	5.31
10	1.42	2.45	1.91	6.75	6.83	5.93
15	1.73	1.27	1.65	6.02	6.57	6.78
20	2.37	1.69	2.31	5.02	6.60	6.51
25	1.27	0.92	1.92	3.57	4.23	2.47
30	0.98	0.95	1.24	4.86	6.47	5.26
35	2.75	1.48	2.44	5.44	5.27	6.50
40	2.58	0.69	1.25	5.05	6.85	7.24
50	1.16	2.29	2.42	6.57	7.17	7.75

Figure 16.

Variation of the accelerometer scale factor with temperature

Figure 17.

Variation of the accelerometer biases with temperature

Figure 18.

Variation of the gyroscope scale factor with temperature

Figure 19.

Variation of the gyroscope biases with temperature

To evaluate the performance of the proposed calibration method, we compare it with the six-position method used in [5] which is the most commonly used. This requires the IMU to be mounted on a levelled table with each sensitive axis pointing alternately up and down. For a triad of orthogonal sensors, this results in a total of six positions. The thermal calibration is implemented by the traditional method for the same temperature. Next, the calibrated measurements of the static IMU are compared by two calibration methods. The results are shown in Figures 20–21. By comparison with the traditional method, the proposed method can improve the measurement precision by an order of magnitude.

Figure 20.

The measured accelerations with zero-inputs

Figure 21.

The measured angular velocities with zero-inputs

5.3 Discussion

Light-weight and low-cost MIMUs are widely used. However, these units need calibration for an accurate measurement solution. The computed biases, scale factors, package misalignment angles and sensor-to-sensor misalignment angles are estimated in a stable temperature through the optimal calibration method. By observing the estimated parameters of the MIMU, we can find out the differences between sets of sensors, even though they belong to the same type. In fact, the parameters' values are not identical to the specifications' values. Meanwhile, it is evident that the accelerometer triad contains smaller orthogonalization errors than the gyros, since the accelerometer triad consists of one three-axis sensor while the gyroscope triad has three single-axis sensors, which are not easy to install orthogonally.

We get the misalignment angles in the optimal calibration at a stable temperature, and then estimate the biases and scale factors of the MIMU at different temperatures by thermal calibration. The experiments have shown that the random errors' variances of the measurements are different. It is useful for the estimation of the biases and scale factors by the Kalman filter. As observed in Figure 17 and Figure 19, there is an almost linear relationship between the biases and temperatures. Hence, a linear interpolation formula can be deduced for calibration in the temperature range. As observed from Figure 16 and Figure 20, the scale factors have an ambiguous relationship with temperature, but the interpolation polynomial can be obtained. The temperature drift errors can be reduced by thermal calibration.

6. Conclusion

The calibration of MIMU is an important phase in improving performance. This paper has presented an optimal calibration method for the MIMU, together with deduction, simulation and experimental results. We analyse the random errors by overlapping Allan variance, and the white noise is taken account of by the calibration. The general calibration model includes such calibration parameters as scale factors, biases, package misalignment error angles and sensor-to-sensor misalignment angles. A new concept for a calibration method is proposed. Moreover, the new calibration scheme is presented by the D-optimal method, the new calibration algorithms are designed with Kalman filters, which are used to obtain optimal estimates of the calibration parameters. Simulations have verified their feasibility and effectiveness. Thermal calibration has been accomplished because the biases of accelerometers and gyroscopes vary significantly with temperature. The experiment results have demonstrated that the new approach outperforms traditional methods. In future, the calibrated MIMU will be used for remotely-operated vehicle navigation for a nuclear power plant.