IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 61, NO.

9, SEPTEMBER 2012 2501

Three-Axial Accelerometer Calibration Using

Kalman Filter Covariance Matrix for Online
Estimation of Optimal Sensor Orientation
Tadej Beravs, Janez Podobnik, and Marko Munih, Member, IEEE

Abstract—Inexpensive inertial/magnetic measurement units can [4], and measurement of human body kinematics [5], where
be found in numerous applications and are typically used to IMUs can replace optical measurement systems, which mea-
determine orientation. Due to the presence of nonidealities in sure body part orientations. [6]. Consisting of accelerometers,
measurement systems, the calibration of the sensor is thus needed
to determine sensor parameters such as bias, misalignment, and gyroscopes, and (optionally) magnetometers, IMUs can achieve
gain/sensitivity. In this paper, an online automatic calibration good dynamic specifications with a relatively low investment.
method for a three-axial accelerometer is presented. Parameters Several IMUs are commercially available; however, custom
are estimated using an unscented Kalman filter. The sensor is developed IMUs have some advantages such us small size,
placed in a number of different orientations using a robotic arm. which allows integration in various applications, custom wire-
These orientations are calculated online from the parameter co-
variance matrix and represent estimated optimal sensor orienta- less connectivity, and open architecture, which allows different
tions for parameter estimation. Numerous simulations are run to modifications and implementations of algorithms. However,
evaluate the proposed calibration method, and a comparison is similar to any measurement systems, IMUs also suffer from
made with an offline least mean squares calibration method. The numerous disadvantages such as sensor misalignment, large
simulation results show that calibration with the proposed method offset, nonlinearity, drift, and random noise.
results in higher accuracy of parameter estimation when using
less than 100 iterations. The proposed calibration method is also These disadvantages are generally addressed using sensor
applied to a real accelerometer using a low number of iterations. calibration. Several offline calibration methods exist. One sim-
The results show only slight (less than 0.4%) changes in parame- ple method proposes the calibration of two main parameters
ter values between different calibration runs. The proposed cali- with manual sensor movement in six different orientations with
bration method provides an accurate parameter estimation using a relatively simple mathematical algorithm (the sum of output
a small number of iterations without the need for manually pre-
defining orientations of the sensor, and the method can be used in signals is equal to the gravity vector) [7]. Similar approaches
combination with other offline calibration methods to achieve even that demand several different sensor orientations are described
higher accuracy. in [8], where all three parameters are determined through the
Index Terms—Accelerometer calibration, orientation determi- Levenberg–Marquardt algorithm, similar to [9], where the pa-
nation, sensor parameter estimation, unscented Kalman filtering rameters are determined using the Newton iterative arithmetic.
(UKF). More sophisticated methods are described in [10], where sensor
parameters are estimated using optical alignment and a least
I. I NTRODUCTION mean squares algorithm, and in [11], the reference orientation
of the sensor is also included in determining the sensor param-

M ICROELECTROMECHANICAL inertial measurement

units (IMUs) are inexpensive lightweight sensors that
are used for orientation estimation in numerous applications.
eters. One method where the robot arm is used to position the
sensor to a known predefined orientation is presented in [12],
where the parameters of the sensor (including the alignment
They can be found in inertial navigation systems [1], [2], angle of the robot in the gravitational field and the alignment
robotics [3], the automotive industry, analysis of daily activities between the sensor and robot end effector) are again determined
using least mean squares methods.
However, because several of the factors that contribute to
sensor errors are time varying (e.g., temperature), initial offline
calibration cannot completely negate their effects. Thus, an
Manuscript received September 20, 2011; revised December 9, 2011; ac- online calibration procedure could potentially achieve higher
cepted December 11, 2011. Date of publication March 5, 2012; date of accuracy. Compared with offline calibration, where parameters
current version August 10, 2012. This work was supported in part by the
European Union through the EVRYON Collaborative Project STREP un- are estimated using different mathematical algorithms after ob-
der Grant FP7-ICT-2007-3-231451 and by the Slovenian Research Agency. taining all measurements from the sensor, the online calibration
The Associate Editor coordinating the review process for this paper was estimates parameters during each iteration, and after the last
Dr. Subhas Mukhopadhyay.
The authors are with the Laboratory of Robotics, Department of Measure- iteration, the estimation of parameters is completed. Online
ment and Robotics, Faculty of Electrical Engineering, University of Ljubljana, approaches have been demonstrated in [13] and [14] using
1000 Ljubljana, Slovenia. sophisticated hardware, but the different sensor orientations
Color versions of one or more of the figures in this paper are available online
at [Link] needed for the calibration must be predefined in both cases.
Digital Object Identifier 10.1109/TIM.2012.2187360 An appropriate orientation must be chosen, or a large number

0018-9456/$31.00 © 2012 IEEE

2502 IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 61, NO. 9, SEPTEMBER 2012

of random orientations must be determined to accomplish an

accurate calibration.
In this paper, we present an automatic calibration method of
the accelerometer, where parameters and orientations are esti-
mated by an unscented Kalman filter (UKF), and a robotic arm
is used to place the sensor in the calculated orientation. Unlike
[12], where the sensor is placed in a large number of manually
predefined orientations using a robotic arm and the parameters
are calculated offline, our proposed method uses online param-
eter estimation without the need for predefined orientations,
because they are calculated and used during calibration. The
method described in [13] uses online parameter estimation;
however, the orientations of the sensor must be predefined.
Fig. 1. IMU that consists of a three-axis gyroscope, a three-axis magnetome-
The method is used to determine all three main parameters ter, and a three-axis accelerometer and a wireless module with a dual-chip
(gain, misalignment, and bias) together with the alignment antenna. The size of the IMU is 30 × 20 × 5 mm without a battery.
angles of the robot in the gravitational field and alignment
angles between the sensor and robot end effector (because electronics three-axis accelerometer; and 3) a Honeywell three-
the flange and sensor board are not perfectly aligned). The axis magnetometer. The gyroscope has selectable full-scale
proposed method repeatedly uses the covariance matrix decom- ranges of ±250◦ /s, ±500◦ /s, ±1000◦ /s, and ±2000◦ /s and
position for estimation of maximal sensitivity axis (CEMS) to software-selectable low-pass filters. Each axis is represented
estimate the next orientation in which the sensor should be with 16 b. The gyroscope also measures temperature for ad-
placed for optimal parameter estimation. This condition causes ditional software compensation. The sampling rate of the gyro-
fast method convergence. The sensor is thus placed in a small scope is 1000 Hz. Similar to the gyroscope, the accelerometer
number of automatically determined orientations, eliminating also offers a selectable range of ±2, ±4, and ±8 g and has 16-b
the need for a large number of predefined orientations and, this output per axis. It offers software selection of high-pass filters
way, allowing faster calibration compared to methods where and sampling rates. The highest possible sampling rate of the
sensor orientations must be predefined and the manipulation of accelerometer is 1000 Hz. The magnetometer has a selectable
the sensor is manually done. Because the sensor is placed in range of ±0.88–±8.1 G. It uses an internal 12-b analog-to-
orientations that allow the most effective parameter estimation digital converter and has a significantly lower sampling rate
and all the data can be recorded, offline methods can also be compared with the other two sensors. The sampling rate of
applied later for parameter estimation. the magnetometer can be selected from 0.75 Hz to 160 Hz.
The CEMS calibration approach can be applied for the ac- Thus, the maximum sampling rate of the IMU system is 160 Hz
celerometer or the magnetometer. The only difference between when data from all three sensors, including the magnetometer,
the two sensors is in the initial description of the gravita- are simultaneously acquired. All sensors are connected to an
tional and magnetic fields. However, the magnetic field is very interintegrated circuit (I2C) bus with a maximum data transfer
sensitive to environmental noise, and a homogenous magnetic rate of 222 kb/s. Each sensor provides 6 B of information
field is needed for successful calibration. The CEMS calibra- (2 B per axis), for a total of 18 B. The theoretically attainable
tion method is thus applied here only to the accelerometer, data transfer rate of the I2C communication protocol is 1.2 kHz,
because the magnetic field that surrounds the robot arm is not but the maximum data transfer rate is set to 300 Hz due to
homogenous. limitations of the wireless transceiver module that provides the
This paper is organized as follows. The developed wireless connection to the central unit. The IMU itself (without battery)
IMU system and the corresponding mathematical model of the has a size of 30 × 20 × 5 mm and is shown in Fig. 1. The
sensor system in conjunction with the robot arm are described battery is placed away from the IMU to avoid interference with
in the first part of Section II. Parameter estimation with UKF is the magnetometer.
described together with the method for determining the sensor 2) Data Transmission and Central Unit: The IMU is con-
orientation using singular value decomposition (SVD) in the nected to a central unit through a 2.4-GHz wireless transmission
second part of Section II. The simulation and measurement pro- system. On the IMU side, a ZigBit wireless transceiver is used.
cedures are described at the end of the section. Simulation and On the receiving side, a powerful Atmel ZigBit receiver is used,
measurement results are presented in Section III, and a detailed because there are no constraints on power consumption and
discussion is given in Section IV. Section V summarizes the size. This receiver has an amplified port for an external antenna
proposed calibration method and the contributions of this paper. and allows a working range of more than 15 m. The receiver is
connected through the serial peripheral interface bus (SPI) to a
central unit, which can simultaneously receive data from up to
II. M ETHODS 10 IMUs at a frequency of 300 Hz and transfer it to a personal
computer through the User Datagram Protocol (UDP). Each
A. Hardware Design
data package is equipped with the time stamp that is generated
1) IMU: The IMU consists of the following three digital on the IMU side together with the checksum. The data from the
sensors: 1) an Invensense three-axis gyroscope; 2) an STmicro- sensor are written as 16-b unsigned integers and are added to
BERAVS et al.: ACCELEROMETER CALIBRATION USING KF MATRIX FOR ESTIMATION OF SENSOR ORIENTATION 2503

the data package. The checksum is then verified by the central

unit, whereas the sensor data are transformed to real numbers
on a personal computer.

B. Kinematic Model of the Sensor and Robotic Arm

A basic mathematical model of a three-axis accelerometer
that includes scaling, misalignment, and bias parameters can be
described as

y =s·T·u+b+N (1)

where vector v represents the output of the sensor for the x-, y-,
and z-axes, vector s = [sx sy sz ] denotes the sensitivity factor
for each axis, and matrix T is described as
⎡ ⎤
1 0 0
T = ⎣ cos α 1 0 ⎦ (2)
cos β cos γ 1

where α, β, and γ represent misalignment angles, vector b =

[bx by bz ] denotes the bias, N represents the noise, and vector
u = [ux uy uz ] denotes the gravitational-field projection on
sensor axes [15]. Because the accelerometer is stationary during
calibration, the only acceleration measured by the sensor is
due to the gravity. The sensor is therefore calibrated in the
range of ±1 g and not by its full-scale range; however, this
Fig. 2. Complete transformation of the gravity vector. Re_b presents the
condition does not represent an issue, because the sensor is transformation between the gravitational field and robot base, R6 presents the
used to determine the orientation of the IMU relative to the transformation between the robot base and the robot end effector, and R6_i
gravitational field. This mathematical model is only a rough presents the transformation between the robot end effector and the IMU.
estimate of a real accelerometer model, because nonlinearity,
temperature drift, and other nonidealities are not considered.
A precise orientation of the sensor can be determined when are determined as follows:
the accelerometer is attached to the robotic arm. A trans- ⎡ ⎤
1 0 0
formation matrix R6 from the robot base frame, denoted as
RotX = ⎣ 0 cos ϕx − sin ϕx ⎦ (5)
coordinate system Ob in Fig. 2, to the end effector OR6 can
0 sin ϕx cos ϕx
be calculated from the robot joint angles using the Denavit–
Hartenberg table. A detailed description of the procedure can ⎡ ⎤
cos ϕz − sin ϕz 0
be found in [12]. Assuming that the robot is perfectly leveled RotZ = ⎣ sin ϕz cos ϕz 0 ⎦. (6)
with the gravitational field, an ideal transformation between the 0 0 1
gravitational field and the projection of the gravitational field
on the sensor uideal can be calculated by Similar to the transformation between the gravitational field
and the robot base, a transformation between the robot end
m = R6 · uideal (3) effector and accelerometer must also be taken into account due
to possible installation errors, because the accelerometer sensor
where m = [1 0 0] represents the unit vector of gravity. How- is not perfectly aligned with the circuit board, and the circuit
ever, perfect alignment of the robot base frame in the grav- board is not perfectly aligned with the end effector. Thus, the
itational field is difficult to achieve. Thus, a transformation transformation matrix R6_i , where φx and φz denote rotation
matrix between the gravitational field, denoted as coordinate over the x- and z-axes, can be described as
system Oe , and the robot base, denoted as coordinate system
Ob , must be taken into account. A transformation matrix can be R6_i = RotZ(φz ) · RotX(φx ). (7)
described as
With both rotational matrices known, a transformation be-
Re_b = RotZ(ϕz ) · RotX(ϕx ) (4) tween the gravitational field and the real projection of the
gravitational field on the sensor u can be calculated by
where ϕx and ϕz denote rotation angles around the x- and
z-axes in coordinate system Oe . Functions RotZ and RotX m = Re_b · R6 · R6_i · u. (8)
2504 IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 61, NO. 9, SEPTEMBER 2012

A complete transformation of the gravity vector is presented of the distribution of ŵk̄ and is usually set to 2 for Gaussian
in Fig. 2. With the transformation matrix specified, the output (c) (m)
distributions. The weights wi and wi are calculated using
of the accelerometer can be described as
(m) λ
w0 =
y = s · T · R−1 −1 −1
6_i · R6 Re_b · m + b + N. (9) L+λ
(c) λ
w0 = + 1 − αkf
2
+ βkf
L+λ
1
C. Parameter Estimation wi (c) =wi (m) = . (15)
2(L + λ)
The UKF is an extension of the traditional Kalman filter
for the estimation of nonlinear systems that attempt to re- The matrix χk|k−1 can be described as
⎡ ⎤
move some of the shortcomings of the extended Kalman filter s0 s1 ··· s2L
(EKF) in the estimation of nonlinear systems. For parameter ⎢ t0 t1 ··· t2L ⎥
⎢ ⎥
estimation, the EKF can be used, because the computation ⎢ b0 b1 ··· b2L ⎥
χk|k−1 = ⎢ (6_i) ⎥ (16)
time of the UKF is greater than the computation time of the ⎢r 0 r(6_i) 1 · · · r(6_i) 2L ⎥
⎣ (e_i) ⎦
EKF. However, because there are no limitations with regard r 0 r(e_b) 1 · · · r(e_b) 2L
to computation time and it has been shown that the UKF n0 n1 ··· n2L
outperforms the EKF in numerous examples, the UKF was
chosen for parameter estimation. More detailed discussion of where vector s0 = ŵk̄,(1...3) consists of the first three ele-
the UKF can be found in [16]–[18]. The UKF uses deterministic ments of vector ŵk̄ . Vectors si = s0 + γ Ŝw̄k i and sL+i =
sampling to approximate the state distribution. The unscented s0 − γ Ŝw¯k i , where i = 1 . . . L are calculated by adding the
transformation uses a set of sample or sigma points that are sigma-point value that was calculated from the ith column of
determined from the a priori mean and covariance of the state. the covariance matrix. A similar approach is applied to vectors
The sigma points are propagated through the nonlinear system. (6_i) (6e_b)
ti , bi , ri and ri . Noise vectors are defined as n0 =
The posterior mean and covariance are then calculated [0 0 0], ni = +γ Ŝw̄k i and nL+i = −γ Ŝw̄k i , where i = 1 . . . L.
from the propagated sigma points. Parameter estimation equa- The output of the sensor model is described as
tions for the UKF are similar to the state estimation. This
section expounds on the differences. yi = si · Ti · R−1 −1 −1
6_i i · R6 k · Re_b i · m + bi + ni (17)
The filter is initialized with the predicted mean and covari-
ance of the parameters, i.e., where values for the matrix Ti are derived from the vector
ti , values for the matrix R6_ii are derived from the vector
ŵ0 = E{w}
  (10) r(6_i) i , and values for the matrix Re_bi are derived from the
Pŵ0 = E w − ŵ0 )(w − ŵ0 )T (11) vector r(6e_b) i . Values for the matrix R6k are obtained from
the orientation of the robotic arm. The expected measurement
where E{ } is the expectation operator, (w − ŵ0 ) is the esti-
values are determined in the matrix Y k|k−1 as
mation error of initial value, w is the unknown true parameter,
and ŵ0 is the estimated initial parameter value. The UKF time Y k|k−1 = [ y0 · · · y2L ]. (18)
update is described as
The measurement mean d̂k̄ and the measurement covariance
ŵk̄= ŵk−1 (12)
Pd̃k are calculated based on the statistics of the expected
P̂w¯k = ηn Pwk−1 + Rwk (13) measurements as
where parameter vector ŵk̄ = [sx sy sz α β γ bx by 2L
bz ϕx ϕz φx φz ] is updated using previous values, and the d̂k̄ = wi (m) Y i,k|k−1 (19)
covariance matrix P̂w¯k is calculated by scaling the previous i=0
value with ηn ∈ (0, 1] and by adding fixed system process noise 2L   T
Rwk . The sigma points χk are calculated from the values of the Pd̃k = wi (c) Y i,k|k−1 − d̂k Y i,k|k−1 − d̂k̄ + Rek .
i=0
mean and covariance of the parameters, i.e.,
 (20)
χk|k−1 = ŵk̄ ŵk̄ + γ Ŝw¯k ŵk̄ − γ Ŝw¯k (14)
The cross-correlation covariance Pwk dk is calculated using

where Ŝw¯k = P̂w¯k is a square root of the covariance matrix

2L
(c)   T
Pwk dk = wi χi,k|k−1 − ŵk̄ Y i,k|k−1 − d̂k̄ + Rek .
T
of wk , P̂w¯k such that
√ P̂w¯k = Ŝw¯k Ŝw¯k2 . Scaling parameters
i=0
are defined as γ = L + λ and λ = αkf (L + κ) − L, where (21)
L denotes the state dimension. The constant αkf determines The Kalman gain matrix is the product of the cross-correlation
the spread of the sigma points around ŵk̄ and is usually set to and measurement covariances, i.e.,
1e − 4 ≤ αkf ≤ 1. κ is a secondary scaling parameter and is
usually set to 0, and βkf is used to incorporate prior knowledge Kk = Pwk dk P−1
d̃
. (22)
k
BERAVS et al.: ACCELEROMETER CALIBRATION USING KF MATRIX FOR ESTIMATION OF SENSOR ORIENTATION 2505

The measurement update equations are given as follows:

w̃k = ŵk̄ + Kk (dk − d̂k̄ ) (23)

Pwk = Pw̄k − Kk Pd̃k KT k (24)

where dk is the measurement from the real sensor or a

simulated output of the sensor, where predefined parameters
are used.

D. Determination of Sensor Orientation

During the parameter estimation, the sensor must be placed
in different orientations to acquire an appropriate set of mea-
surements for successful parameter estimation. In the proposed
CEMS algorithm, the orientation is chosen to position the sen-
sor in orientation, in which the maximal sensitivity is achieved
for parameters with the largest variance. This orientation can
be determined from the covariance matrix Pwk . The Kalman
filter returns the estimation of the posterior mean state ŵk
and error covariance Pwk . The posteriori error covariance
Pwk is segmented into two covariance submatrices that rep-
resent the posteriori error covariance matrices of bias and gain. Fig. 3. Simplified two-degree-of-freedom example of the error ellipse of the
covariance matrix of parameter estimation error. (a) Initial error ellipse. (b) and
The posteriori error covariance matrices of bias and gain are (c) Intermediate steps. (d) Final error ellipse. Axes s1 and s2 are the principle
used as covariance matrices of the parameter estimation error, axes of error ellipse, where principle axis s2 has a smaller variance.
defined as
  [σ1 σ2 σ3 ]. Singular values are associated with the variance.
Pwk = E (wk − ŵk )(wk − ŵk )T (25) Singular value σ3 is associated with the lowest variance, and
thus, a unit vector u3 corresponds to the principle axis with the
where E{ } is the expectation operator, (wk − ŵk ) is the lowest variance of the covariance of parameter estimation error.
estimation error, wk is an unknown true parameter value, and An intuitive interpretation of the proposed SVD approach is
ŵk is an estimated parameter value. given by the principal component analysis (PCA). PCA uses
The covariance matrix of the parameter estimation error is an orthonormal transformation to transform the original space
positive semidefinite and a symmetric matrix and can therefore into a new one, where the first axis points in the direction
be diagonalized using an orthonormal basis. The unit vectors of the maximum variance and the subsequent axes are ranked
of the orthonormal basis used to rotate the covariance matrix according to the variance with the final axis pointing in the
are the eigenvectors of the covariance matrix and form the direction of the lowest variance [20].
principle axes of an error ellipse. The values of the diagonalized This methodology is used with a stationary accelerometer.
covariance matrix are the eigenvalues of the covariance matrix The measurement thus corresponds only to the projection of
and correspond to the variances of the decoupled noise con- the gravitational field. The orientation vector of the sensor is
tributions in the direction of the corresponding principle axes calculated from the singular vector u3 . The singular vector
of the error ellipse. Fig. 3 shows the simplified two-degree-of- u3 corresponds to the principle axis with the lowest variance
freedom example of error ellipse with two principle axes s1 and expressed in the coordinate system of the IMU. Because the
s2 . Fig. 3(a) shows the initial error ellipse, where a large initial axis is needed to control the robot, the principle axis must be
value of variance is chosen, and both principle axes have same transformed into the coordinate system of the endpoint of the
variance values. Fig. 3(b) and (c) shows the intermediate steps. robot with transformation, i.e.,
Fig. 3(d) shows the final error ellipse, where the parameter
estimation error is minimized, and variances are approximately ue = R6_i · u3 (27)
equal for both ellipse principle axes s1 and s2 .
SVD can be used to decompose the covariance matrix into an where ue denotes the orientation vector on the robot end
orthonormal basis and a diagonal matrix. The SVD algorithm is effector.
applied to each of the covariance matrices of parameter estima- After applying the orientation with robot motion, the sensor
tion errors [19]. Because covariance matrix Pw is a positive- is aligned in orientation, which will maximize the sensitivity for
semidefinite symmetric matrix, the following decomposition is the sensor axis with the largest variance, and the sensitivity will
obtained for a selected parameter: be lowest for the axis with the lowest variance. The principal
axis with the lowest variance is positioned to be perpendicular
Svd(Pwpar ) = U · Σ · UT . (26) to the gravitational field so that the performing rotation around
this axis will align the other two principle axes with the
U = [u1 u2 u3 ] is an orthonormal basis matrix of singular gravitational field. The initial orientation is set to the value for
vectors, and matrix Σ is a diagonal matrix of singular values which the principle axis with the largest variance is aligned with
2506 IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 61, NO. 9, SEPTEMBER 2012

Fig. 4. Initial sensor orientation is noted with axes xe , ye , and ze . After one
iteration is completed, the sensor must reach the new orientation noted with
axes xe , ye , and ze . Additional rotation is applied around the ze -axis.

the gravitational field. In Fig. 3(b), the robot will position the
sensor in the orientation for which the principle axes s1 of the
Fig. 5. Flowchart of the simulation procedure of the calibration method.
error ellipse will be aligned with the gravitational field. After
several steps [see Fig. 3(c)], the variance will decrease, and
principle axes of error ellipse will move to different orienta- E. Simulation and Measurements
tions, and therefore, the robot will reorient the sensor to align Simulation is used to verify the kinematic model and pro-
the principle axes s1 with the gravitational field. The final result posed procedure, because the true parameters of the sensor are
is a sequence of movements of rotations around the principle not known. All sensor parameters are manually predefined in
axis, which maximizes the sensitivity of the sensor axis with the the kinematic model and later compared with the simulation
largest variances of the parameter estimation error. After each results, thus allowing us to validate the calibration method. The
new measurement, a new principle axis is computed, and the simulation is built and run in MATLAB. Because the calibration
movement of the robot is updated with the new axis of rotation method is also based on the movement of the robot arm, a
to reduce the variance along the axis with the largest uncertainty simulation of the robot must also be included. The simulation
(see Fig. 4). process can be segmented into the following three parts.
Singular values σi are used to determine the validity of the
• The output of the sensor is calculated (simulated) using
estimated parameters gathered from the UKF filter, because
manually predefined parameters and known rotational ma-
they represent the dispersion around the associated axis. A
trices, i.e., R6 , Re_b and R6_i .
validation criterion for the selected parameter estimation is
• The calculated output of the sensor is fed into the UKF
presented by
algorithm. The orientation result given by the UKF kine-
3 · σ3 matics is described as a unit vector. Estimated parameters
Cpar = . (28)
3 are temporarily stored and used for the next iteration.
σi
i=1
• A fixed rotation is applied over a unit vector that results
in a 3 × 3 rotational matrix that represents the robot end
Three criterion functions Cpar are calculated, because three effector R6 . With the orientation matrix known, the next
parameters are determined through calibration. The closer Cpar simulation step can be performed.
is to 1, the lower the largest variance is compared to the sum of
Because initial values are needed by the UKF, they are set
variances. A value of 1 also implies that the variance is lowest
close to ideal values with a small offset. For example, the initial
as possible, because there is no axis that would further reduce
parameter for bias is set as ba = [0.15 0.2 −0.12]. Although
the variance. This case is shown in Fig. 3(d), where variances of
the ideal parameter is ba = [0 0 0], a small offset allows the
both principle axes of error ellipse are approximately equal. The
parameters to more quickly converge to the true value. After
criterion functions are also used as a weight for determining the
numerous runs of the simulation, the UKF parameters are
rotation axis of the sensor by
adjusted to ensure rapid convergence of the criterion function.
u = (1 − Cb ) · ue_b + (1 − Cs ) · ue_s (29) The flowchart of the calibration is shown in Fig. 5.
Once the simulation parameters are determined, 300 simula-
where u represents the axis of sensor rotation, Cb and Cs tion runs are performed. Because the model of the sensor and
represent the criterion functions of bias and gain/sensitivity, the UKF algorithm have a fixed noise parameter, different simu-
and ue_b and ue_s represent the estimated axis of rotation for lation runs output different sensor and parameter estimates. The
both parameters. When all criterion functions are close to 1, the dispersion of the parameter values around the true predefined
calibration procedure can be completed, because the variance value can be used to evaluate the CEMS calibration method.
is the lowest, and further measurement will not improve the After the simulations are successfully completed, the pro-
estimation of parameters. posed calibration method is applied using the IMU described
BERAVS et al.: ACCELEROMETER CALIBRATION USING KF MATRIX FOR ESTIMATION OF SENSOR ORIENTATION 2507

TABLE I TABLE II
S IMULATION R ESULTS OF E STIMATING G AIN , M ISALIGNMENT, AND B IAS S IMULATION R ESULTS OF E STIMATING A NGLES IN C OORDINATE
PARAMETERS W ITH THE P ROPOSED C ALIBRATION M ETHOD U SING 400 S YSTEMS Oe and OR6 W ITH THE P ROPOSED C ALIBRATION M ETHOD
I TERATIONS . T HE F IRST C OLUMN P RESENTS THE P REDEFINED VALUES , U SING 400 I TERATIONS . T HE F IRST C OLUMN P RESENTS THE
THE S ECOND C OLUMN P RESENTS THE C ALCULATED M EAN VALUES , THE P REDEFINED VALUES , THE S ECOND C OLUMN P RESENTS THE
T HIRD AND F OURTH C OLUMNS P RESENT THE M INIMUM AND M AXIMUM C ALCULATED M EAN VALUES , THE T HIRD AND F OURTH C OLUMNS
VALUES , THE F IFTH C OLUMN P RESENTS THE M EDIAN VALUES , AND P RESENT THE M INIMUM AND M AXIMUM VALUES , THE F IFTH C OLUMN
THE S IXTH C OLUMN P RESENTS THE S TANDARD D EVIATION P RESENTS THE M EDIAN VALUES , AND THE S IXTH C OLUMN
P RESENTS THE S TANDARD D EVIATION

in Section II and a six-axis Epson PS3 robot. The IMU is

tightly attached to the aluminum flange that is bolted to the
robot end effector. Data from the sensor are wirelessly trans-
mitted to the receiver board, which is connected to a personal
computer through the UDP. Data from the IMU are acquired
and transferred to the UKF using MATLAB/Simulink. Once the
UKF calculation is done, the new orientation of the sensor must
be transmitted to the robot. Because the Epson robot accepts
orientation in values of angles over the x-, y-, and z-axes, the
Fig. 6. Histogram of the difference of estimated gain from the known true
orientation matrix must be transformed into these three angles. gain value for the x-axis.
Because position is not relevant for accelerometer calibration,
the position can be changed to achieve the desired orientation. TABLE III
This approach cannot be done for magnetometer calibration, S TANDARD D EVIATIONS FOR G AIN , M ISALIGNMENT, AND B IAS
PARAMETERS U SING 500 S IMULATION RUNS . VALUES A RE C ALCULATED
because it is difficult to ensure a constant magnetic field in U SING THE D IFFERENCE B ETWEEN THE E STIMATED AND K NOWN T RUE
the surroundings of the robot. The three orientation angles VALUES , W HICH A RE R ANDOMLY C HANGED B ETWEEN
are received by the Epson robot through the Transmission D IFFERENT S IMULATION RUNS
Control Protocol/Internet Protocol (TCP/IP). Once all data are
received, the robot moves to the specified orientation with a
low speed. This case avoids any vibrations that could occur
during movement, because the robot arm is not perfectly rigid.
After the robot reaches the desired orientation, a signal flag that
indicates that the robot is stationary is sent to MATLAB, and
the new acquisition of the sensor data can commence.
Similar to the simulation, multiple measurements/
calibrations are performed with the same IMU to evaluate the
calibration method by comparing measured sensor parameters
of the sensor. A fixed number of iterations (400) are used Similar to Table I, Table II presents the estimated values of
for each measurement. This number is determined from the angles in coordinate systems Oe and system OR6 . The first
parameter criterion function during simulation. column presents predefined values of angles, the second column
presents mean values, the third and fourth columns present
minimum and maximum values, the fifth column presents the
III. R ESULTS median, and the sixth column presents the standard deviation.
Standard deviations of parameter estimations are determined
A. Simulation
by 500 simulation runs, and the values of predefined parameters
Evaluation of the method is done by running 100 simulations. are randomly changed for each simulation. The gain parameters
Predefined parameters of the sensors are listed in Table I, first are in the range from 0.9000 to 1.1000, the misalignment
column, whereas the second column presents the mean values parameters are in the range from 1.4708 to 1.6708, and the bias
of calculated parameters within all simulations, the third and parameters are in the range from −0.1500 to 0.1500. The data
fourth columns present the minimum and maximum values of obtained are used for the calculation of the standard deviation of
parameters that occurred during evaluation, the fifth column parameters for each axis. The differences between the estimated
presents the median value, and the sixth column presents the gain from the known true gain value for the x-axis are presented
standard deviation. The results presented in Table I are mea- in the histogram in Fig. 6. Standard deviations of parameters
sured with 400 iterations. gain, misalignment, and bias are presented in Table III.
2508 IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 61, NO. 9, SEPTEMBER 2012

Fig. 7. Mean and maximum errors of gain and misalignment parameter esti- Fig. 9. Mean and maximum gain and misalignment errors using the least
mations using the proposed calibration method. Dashed lines present maximum mean squares method. Dashed lines present maximum errors that occurred
errors that occurred during simulations, and solid lines present mean relative during simulations, and solid lines present mean relative errors.
errors.

Fig. 10. Mean and maximum bias offsets using the least mean squares
Fig. 8. Mean and maximum bias offsets using the proposed calibration method. The dashed line presents the maximum error that occurred during
method. The dashed line presents the maximum error that occurred during simulations, and the solid line presents the mean relative error.
simulations, and the solid line presents the mean relative error.
TABLE IV
To present an overview of how the number of measurements E STIMATED PARAMETERS OF THE R EAL IMU U SING F IVE D IFFERENT
M EASUREMENTS W ITH 400 I TERATIONS
or iterations influences the accuracy of the parameter estima-
tion, a new series of simulations is run with a variable number of
iterations within the range of 20–500 iterations. Fig. 7 presents
the mean and maximum relative error of gain and misalignment
parameters as solid and dashed lines. The error is calculated by
running 100 simulations for each selected number of iterations.
The bias, misalignment, and gain parameters are calculated
for each axis, and the success of the calibration is determined
by the worst parameter estimated. Thus, only the maximum
relative error that occurred on any of three axes during a single
simulation run is used for the calculation of a mean relative proposed calibration method. The maximum number of offline
error of 100 simulation runs. The maximum relative errors (more than 5000) iterations is used to calculate the values
that occurred during simulation runs at different numbers of of parameters. For each number of movements, 100 offline
iterations are also presented in Fig. 7, dashed lines. iterations are run. The same method for the calculation of the
Because the preset bias parameters are near zero, Fig. 8, mean parameter estimation error is used. The mean parameter
solid line, presents the mean offsets between the calculated and estimation errors and the maximum relative errors that occurred
preset values during 100 simulations at different numbers of during simulation runs at different numbers of iterations are
iterations. Values are calculated using the maximum offset that presented in Fig. 9, dashed and solid lines. The mean and
occurred on any of three axes during simulation runs. Similar to maximum bias offsets at different numbers of iterations are
the previous figure, the maximum values that occurred during presented in Fig. 10.
simulation runs are also presented as a dashed line.
Further evaluation of the CEMS calibration method is done
B. Real IMU
compared with the commonly used least mean squares method,
which is an offline method that requires a different approach. Measurements of a real IMU are performed using a robotic
Because this method cannot set the orientation of the sensor, arm. Because precise parameters of a real sensor are unknown,
a random movement is generated. For better comparison, the Table IV presents the estimated values. Comparison between
number of movements is equal to the number of iterations in our five measurements is made, each with 400 iterations.
BERAVS et al.: ACCELEROMETER CALIBRATION USING KF MATRIX FOR ESTIMATION OF SENSOR ORIENTATION 2509

bias calculated with the proposed calibration method according

to Table I is in the range of 0.0015 g, whereas according to
Table III, the standard deviation can be up to 0.0039 g and is
in the acceptable range. A comparison between the estimated
and real misalignment values cannot be made, because the
manufacturer does not provide this information. However, the
misalignment mean relative error is also lower than 0.5%, and
the standard deviation is up to 0.0136 rad. The minimum and
maximum values noted in Table I, which occurred during 100
different simulation runs, are used to calculate the maximum
gain and misalignment relative error, which is 4.5%, and for the
Fig. 11. Criterion functions of gain, misalignment, and bias parameters during calculation of the maximum bias offset, which is 0.02 g.
400 iterations. The estimation of angle parameters for coordinate system
OR6 have, according to Table II, a mean error of less than
0.02%, and the difference between the minimum and maximum
values does not exceed more than 0.0020 rad. The estimation of
angle parameters for coordinate system Ob have slightly higher
error when taking into account the difference between the
minimum and maximum values, which is up to 0.0140 rad (see
Table II). The rotational matrix R6 is, in simulation, determined
as absolutely accurate; however, when the calibration method
is used on a real robot, the accuracy of this matrix depends on
the robotic arm accuracy, which can be determined from robot
specifications.
Fig. 7 presents the influence of varying numbers of iterations
on parameter estimation accuracy. As expected, the highest
Fig. 12. Values of parameters during 400 iterations. The upper part presents
values of misalignment parameters, the middle part presents the values of gain mean relative error (lower than 1.4%) is achieved with the low-
parameters, and the lower part presents the values of bias parameters. est number of iterations (20). However, the maximum relative
error that occurred during the calibration is up to 10.7%. The
Fig. 11 presents the criterion functions of bias, misalignment, mean offset of the estimated bias is 0.012 g at 20 iterations,
and gain parameter estimations during calibration. The values as shown in Fig. 8. Similar to Fig. 7, the maximum deviation
shown in the figure are calculated from the measurement of a of the bias is much higher than the mean value, i.e., −0.1 g.
real IMU during 400 iterations. When the number of iterations is increased, the gain and
Similar to the Fig. 11, the estimation of the parameter values misalignment mean relative error and bias mean offset slightly
is observed and presented in Fig. 12 during 400 iterations for decrease, whereas the decreases of the maximum relative error
the real IMU calibration. In this figure, the upper part of the and maximum offset are much more notable. The misalignment
plot presents values of misalignment parameters for each sensor and gain maximum relative errors decrease to 5% and 7%,
axis, the middle part presents the values for gain parameters, respectively, whereas the maximum bias offset decreases to
and the lower part presents the values for the bias parameters 0.04 g. Further increase of the iteration number does not
for each sensor axis. The x-, y-, and z-axes are marked with significantly decrease the mean and maximum relative error
solid, dashed, and dotted lines, respectively. or mean and maximum bias offset. There is a convergence of
0.57% for the mean relative error, 5% for the maximum relative
error, 0.0007 g for the mean bias offset, and 0.025 g for the
IV. D ISCUSSION
maximum bias offset.
According to the results presented in Table I, the CEMS cal- The comparison of our method and the offline least mean
ibration method can estimate parameters with a mean relative squares calibration method in Figs. 7 and 9, as well as in Figs. 8
error of 0.5% when 400 iterations are performed. However, and 10, clearly shows that our proposed calibration method
according to Table III, the maximum standard deviation for gain yields much lower errors at a low number of iterations. The
parameter estimation is 0.0096. Because the gain is in a range misalignment and gain maximum relative errors are higher than
of value 1, the relative error of determining the gain parameter 30%, and the mean relative errors are higher than 7% when
is less than 1%. Focusing on the real accelerometer, the gain using the offline least mean squares method. Similarly, the
parameter is within ±10% of the true value according to the maximum bias offset is 0.26 g, and the mean offset is lower
manufacturer’s specification. The accuracy of the estimated than 0.08 g, which is close to the maximum bias offset of our
gain parameter is therefore within the acceptable range, because calibration method. However, the inaccuracy of parameter esti-
it is much higher than the gain accuracy of the uncalibrated mation is due to the low number of random sensor orientations.
accelerometer. Deviation from the ideal value (zero) in the These orientations cannot cover the most influential positions
bias parameter can be within the range of ±0.02 g according where the sensitive axes are aligned with the gravitational field
to the manufacturer’s specification. The mean deviation of the in both directions. Increasing the number of random movements
2510 IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT, VOL. 61, NO. 9, SEPTEMBER 2012

up to 100 greatly reduces the mean and maximum relative parameter estimation accuracy as a function of the number of
errors and offsets. Nonetheless, at 100 iterations, the errors iterations. High accuracy was achieved after a relatively low
of the offline least mean squares calibration method are still number of iterations compared to an offline calibration method
significantly higher than in our proposed calibration method with randomly generated sensor orientations. The proposed
(except for the misalignment maximum relative error, which is calibration method was then applied to a real accelerometer,
lower by 1%). Increasing the number of different orientations where a parameter estimation relative error of less than 0.3%
makes the maximum relative errors converge to 2.4% and 1.1%, was achieved. Although other offline methods could potentially
whereas the maximum offset converges to 0.01 g. The gain and achieve higher accuracies, our approach represents a promising
misalignment mean relative error also decrease by 0.9% and method that can automatically determine appropriate sensor
0.5%, respectively, whereas the mean bias offset is decreased orientations for calibration and thus rapidly produce accu-
to 0.005 g. rate sensor parameters online, without the need for operator
Our proposed calibration method, thus, has an advantage in involvement.
parameter estimation when less than 100 iterations are used,
because the mean and maximum errors are significantly lower
than the errors calculated by the offline calibration method ACKNOWLEDGMENT
with random orientations. Fig. 11 clearly shows that criterion
The authors would like to thank D. Novak for the complete
functions for all three parameters begin to converge to 1 after
overview of this paper and for the grammatical corrections.
50 iterations, which means that the errors are close to their
minimum. In Fig. 12, the parameter values similarly settle
after 50 iterations, and only slight adjustments are made in R EFERENCES
further iterations. Because the optimal sensor orientations are [1] G. Grenon, P. An, S. Smith, and A. Healey, “Enhancement of the inertial
determined by the calibration method, there is no need to navigation system for the Morpheus autonomous underwater vehicles,”
manually find and move the sensor to appropriate orientations, IEEE J. Ocean. Eng., vol. 26, no. 4, pp. 548–560, Oct. 2001.
[2] C. Tan and S. Park, “Design of accelerometer-based inertial navigation
thus automating and shortening the time of calibration. systems,” IEEE Trans. Instrum. Meas., vol. 54, no. 6, pp. 2520–2530,
However, the disadvantage of this method is that it uses Dec. 2005.
relatively expensive equipment for sensor manipulation. Better [3] E. Bachmann, I. Duman, U. Usta, R. McGhee, X. Yun, and M. Zyda,
“Orientation tracking for humans and robots using inertial sensors,” in
accuracy can be achieved with offline calibration methods when Proc. IEEE Int. Symp. CIRA, 1999, pp. 187–194.
a large number of sensor orientations or carefully predefined [4] M. Mathie, B. Celler, N. Lovell, and A. Coster, “Classification of basic
sensor orientations are used. A combination of both methods daily movements using a triaxial accelerometer,” Med. Biol. Eng. Com-
put., vol. 42, no. 5, pp. 679–687, Sep. 2004.
could therefore result in better accuracy of parameter estima- [5] M. Mihelj, “Inverse kinematics of human arm based on multisensor data
tion. However, this approach would extend the total calibration integration,” J. Intell. Robot. Syst., vol. 47, no. 2, pp. 139–153, Oct. 2006.
time, creating a disadvantage compared to our online cali- [6] R. Mayagoitia, A. Nene, and P. Veltink, “Accelerometer and rate gyro-
scope measurement of kinematics: An inexpensive alternative to opti-
bration method, where the parameter estimation is complete cal motion analysis systems,” J. Biomech., vol. 35, no. 4, pp. 537–542,
immediately after the final iteration. Apr. 2002.
The CEMS calibration method, in the future, can also [7] S. Won and F. Golnaraghi, “A triaxial accelerometer calibration method
using a mathematical model,” IEEE Trans. Instrum. Meas., vol. 59, no. 8,
be applied to three-axial magnetometer calibration. Because pp. 2144–2153, Aug. 2010.
the manipulation of the magnetic sensor is not ideal due to [8] F. Camps, S. Harasse, and A. Monin, “Numerical calibration for three-
magnetic-field disturbances around the robotic arm, a modified axis accelerometers and magnetometers,” in Proc. IEEE EIT Conf., 2009,
pp. 217–221.
method that involves a three-axial magnetic coil can be used. [9] J. Wang, Y. Liu, and W. Fan, “Design and calibration of a smart inertial
In this case, the direction of the magnetic field would also be measurement unit for autonomous helicopters using MEMS sensors,” in
determined by the calibration method. Because there would be Proc. IEEE Int. Conf. Mechatron. Autom., 2006, pp. 956–961.
[10] R. Zhu and Z. Zhou, “Calibration of three-dimensional integrated sensors
no need for physically moving the sensor and the change in for improved system accuracy,” Sens. Actuators A: Phys., vol. 127, no. 2,
the magnetic field can instantly be done, this approach could pp. 340–344, Mar. 2006.
represent a very fast method of magnetometer calibration. [11] A. Kim and M. Golnaraghi, “Initial calibration of an inertial measurement
unit using an optical position tracking system,” in Proc. IEEE PLANS,
2004, pp. 96–101.
[12] E. Renk, M. Rizzo, W. Collins, F. Lee, and D. Bernstein, “Calibrating
a triaxial accelerometer-magnetometer-using robotic actuation for sensor
V. C ONCLUSION reorientation during data collection,” IEEE Control Syst. Mag., vol. 25,
no. 6, pp. 86–95, Dec. 2005.
This paper has presented an online automatic calibration [13] P. Batista, C. Silvestre, P. Oliveira, and B. Cardeira, “Accelerometer cal-
method for a three-axial accelerometer. A robotic arm is used ibration and dynamic bias and gravity estimation: Analysis, design, and
experimental evaluation,” IEEE Trans. Control Syst. Technol., vol. 19,
to rapidly place the sensor in a number of different orientations, no. 5, pp. 1128–1137, Sep. 2011.
and the UKF estimates the three main accelerometer parame- [14] H. Kuga, R. da Fonseca Lopes, and W. Einwoegerer, “Experimental static
ters (gain, misalignment, and bias) in each orientation. These calibration of an IMU (inertial measurement unit) based on MEMS,” in
Proc. XIX COBEM, Brasília, DF, Brazil, 2007.
orientations are automatically calculated during calibration us- [15] D. Jurman, M. Jankovec, and R. Kamnik, “Calibration and data fusion
ing the parameter covariance matrix to represent the optimal solution for the miniature attitude and heading reference system,” Sens.
orientations for parameter estimation. Actuators A: Phys., vol. 138, no. 2, pp. 411–420, Aug. 2007.
[16] R. Van Der Merwe, “Sigma-point Kalman filters for probabilistic infer-
Several simulations were performed to evaluate the CEMS ence in dynamic state-space models,” Ph.D. dissertation, Oregon Health
calibration method. Its success was measured by observing Sci. Univ., Portland, OR, 2004.
BERAVS et al.: ACCELEROMETER CALIBRATION USING KF MATRIX FOR ESTIMATION OF SENSOR ORIENTATION 2511

[17] J. Ambadan and Y. Tang, “Sigma-point Kalman filter data assimilation Janez Podobnik received the B.S. degree in elec-
methods for strongly nonlinear systems,” J. Atmos. Sci., vol. 66, no. 2, trical engineering and the Ph.D. degree from the
pp. 261–285, Feb. 2009. University of Ljubljana, Ljubljana, Slovenia, in 2004
[18] M. VanDyke, J. Schwartz, and C. Hall, Unscented Kalman Filtering for and 2009, respectively.
Spacecraft Attitude State and Parameter Estimation, Blacksburg, VA2004. He is currently a Researcher and Teaching Assis-
[19] S. Bonnet, C. Bassompierre, C. Godin, S. Lesecq, and A. Barraud, “Cal- tant with the University of Ljubljana. His research
ibration methods for inertial and magnetic sensors,” Sens. Actuators A: interests include haptic interfaces, real-time control
Phys., vol. 156, no. 2, pp. 302–311, 2009. of robots for virtual-reality-supported rehabilitation,
[20] C. Bishop and S. S. en Ligne, Pattern Recognition and Machine Learning. and sensory fusion techniques.
New York: Springer, 2006.

Marko Munih (M’88) received the Ph.D. degree

in electrical engineering from the University of
Ljubljana, Ljubljana, Slovenia.
From 2004 to 2006, he was the Head of the
Department of Measurement and Robotics, Faculty
of Electrical Engineering, University of Ljubljana,
where he is currently a Full Professor and the Head
of the Laboratory of Robotics. He was a Principal
Investigator for eight European Union (EU) projects.
His early research interests were focused on the
functional electrical stimulation of paraplegic lower
Tadej Beravs received the B.S. degree from the Uni- extremities with surface electrode systems, including measurements, control,
versity of Ljubljana, Ljubljana, Slovenia, in 2010. biomechanics, and electrical circuits. In the last 15 years, his research focus was
He is currently working toward the Ph.D. degree in on robot contact with environment, as well as the construction and use of haptic
the Laboratory of Robotics, Department of Measure- interfaces in the industry and rehabilitation engineering, in combination with
ment and Robotics, Faculty of Electrical Engineer- VR. In the industry, his research interests include specific robot applications in
ing, University of Ljubljana. construction and robots for deburring and measurement tasks, covering laser
He is also a Junior Researcher with the Laboratory technology for the measurement of distance and deviation.
for Robotics, Department of Measurement and Ro- Dr. Munih is the recipient of the Zois Award from the Slovene Ministry of
botics, Faculty of Electrical Engineering, University Science, Education and Sport in 2002 for his outstanding scientific contribu-
of Ljubljana. His research interests include calibra- tions and the Vidmar Award for Best Professor from the Faculty of Electrical
tion methods and development of inertial measure- Engineering, University of Ljubljana, in 2011. He is a corecipient of the 2010
ment units. EUROP/EURON Robotics Technology Transfer Award, Third Prize.