- Email: [email protected]
Methods comparison for attitude determination of a lightweight buoy by raw data of IMU
Accepted Manuscript Methods Comparison for Attitude Determination of a Lightweight Buoy by Raw Data of IMU Arash Joybari, Hadi Amiri, Alireza A. Ardalan, Niloofar K. Zahraee PII: DOI: Reference:
S0263-2241(18)31115-1 https://doi.org/10.1016/j.measurement.2018.11.061 MEASUR 6104
To appear in:
Measurement
Received Date: Revised Date: Accepted Date:
23 May 2018 30 October 2018 18 November 2018
Please cite this article as: A. Joybari, H. Amiri, A.A. Ardalan, N.K. Zahraee, Methods Comparison for Attitude Determination of a Lightweight Buoy by Raw Data of IMU, Measurement (2018), doi: https://doi.org/10.1016/ j.measurement.2018.11.061
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Methods Comparison for Attitude Determination of a Lightweight Buoy by Raw Data of IMU Arash Joybari1, Hadi Amiri2, Alireza A. Ardalan3, Niloofar K. Zahraee4 1
Faculty of Engineering and Sustainable Development, University of Gävle, Sweden; [email protected]. School of Engineering Science, College of Engineering, University of Tehran, Iran; [email protected] 3 School of Surveying and Geomatics Engineering, University of Tehran, Iran; [email protected] 3 School of Surveying and Geomatics Engineering, University of Tehran, Iran; [email protected] 2
Correspondence to: Hadi Amiri ([email protected]) Abstract. Today, one of the most important issues is the determination of instantaneous sea level and distinguishing the Tsunami by floating buoy in the ocean. Usually, gyroscopes are used to measure the angular velocity of a buoy. On the other hand, considering the advancement of various technologies in the field of precise accelerometers, make it possible to use these kinds of sensors for navigation purpose. In this research, stable and optimal methods for determining the orientation of a moving buoy is presented using a combination of the gyroscope, accelerometers, and magnetic sensors data. In order to prove the effectiveness of the proposed methods, the raw data were collected from accelerometers, gyroscopes, and magnetometers of (Xsens MTI-G-700) mounted on a Buoy in coastal waters of Kish Island, Iran. Then, by using the proposed methods, the Euler angles of the buoy are determined, while the Euler angles are derived from the Xsens sensor we are considered as a reference. Based on the results, RMSD for Madgwick algorithm are 0.57°, 0.37°, and 0.50° for Mahony algorithm are 0.56°, 0.37°, and 0.50° and finally for Complementary algorithm is 0.63°, 0.26°, and 2.38° which these values are for roll, pitch, and yaw angles respectively. Thus Mahony algorithm for determining roll and yaw Euler angles is more accurate than other algorithms; however, this differences is negligible compared to the Madgwick algorithm. The Complementary algorithm is less accurate than the other two algorithms, especially for determining the yaw angle of the buoy.
1. Introduction Precise measurement of the angles plays a fundamental role in various applications such as aerospace [1], [2], robotics [3], [4], navigation [5], [6], and human movement analysis [7], [8]. Despite the fact that today there are various technologies for measuring the Euler angles, but the use of inertial systems has become more and more progressive since these systems do not depend on the dynamics of movement and environmental conditions [9]. Indeed, a magnetic and inertial measurement unit (MIMU) consists of gyroscope, accelerometer, and magnetometer which, in the case of an orthogonal installation of sensors along the main axis, has the ability to measure rotational and transient Indeed, an IMU (Inertial Measurement Unit) consists of gyroscope, accelerometer, and magnetometer which, in the case of an orthogonal installation of sensors along the main axis, has the ability to measure rotational and transient movements along the three main axes of the body [5]. A gyroscope is able to measure the sensor’s attitude during the time, if its initial state is known [10], [11]. However, because of 1
the accumulation errors in gyroscope measurements, the calculated Euler angles can have significant errors. Accordingly, accelerometers and magnetometers, which measure gravity and magnetic fields, can be integrated with gyroscopes and thus produce accurate Euler angles. It is worth noting that accelerometers are also exposed to relatively high noise levels, such that sudden motion of the sensor leads to a reduction in the accuracy of the accelerometer's gravity observations. In general, the use of an appropriate integration algorithm for calculating and estimating optimal Euler angles is essential by combining gyroscope, accelerometer, and magnetometer measurements [12]. Typically, the Kalman filter has been used as a wellknown tool in most of the commonly used research in Euler angles determination using raw observations of gyroscope, accelerometer, and magnetometer [13]. The basis of the Kalman process is based on the repetition of linear or non-linear regressions in different epochs. The widespread use of the Kalman filter is evidence of its high precision and performance, although it has some disadvantages. Alternatively, we can use the fuzzy process or fixed filters to determine the Euler angles in order to improve the accuracy of accelerometer observations at low angular velocities and gyroscope observations at high angular velocities [14]. On the other hand, [15] proposed another method based on least squares regression, so that the filter integrates the observations of all inertial sensors in determining the Euler angles. In this research, three methods for integrating the inertial sensor observations are presented in determining the Euler angles while in our previous research just two methods for integrating the inertial sensor observations has been evaluated [16]. The fundamental of this filter is based on the use of quaternion for moving platform angles determination in a three-dimensional space [17]. The Comparison of proposed filters and commercial filters (Xsens Kalman filter) has been investigated. In the next section, the details of the proposed method are discussed. After that, the results of applying the proposed method are presented.
2. Materials and Methods A. Complementary Filter The Kalman filter is very efficient and accurate, but it also has some disadvantages as complexity. The complementary filter is much easier to understand and implementation than the Kalman filter. Most IMUs have 6 degrees of freedom, which means that there are 3 accelerometers and 3 gyroscopes. The IMU has the ability to measure the position and orientation of the platform to which it is attached because an object in space has 6 degrees of freedom. So, if 6 degrees of freedom for an object in space is measured, all problems related to determining of position and orientation of object would be solved. The raw data of the sensors in the IMU unit is not enough to achieve this goal. Observations of accelerometer and gyroscope sensors are used to reach the position of a platform in space. The gyroscope's observations are angular velocities; therefore, by integrating the angular velocity per time unit, the Eulerian angles are obtained. The accelerometer observations calculate the Eulerian angles, though. Accelerometers measure the position and state of the gravity vector (g-force) and calculate the platform orientation using the trigonometric relationship. A major problem with using these two types of observation to reach the platform orientation is that it is not possible to determine the orientation of platform without filtering on them. The problem with accelerometer measurements is that the accelerometers, in addition to gravitational acceleration, measure all the forces on the moving platform. The accelerometer, in addition to measuring gravitational acceleration, also measures noises, so accelerometer observations are reliable only in long-term measurements and require a low pass filter for removing noises. But the gyroscope calculates the Eulerian angles by integrating the angular rate in the unit time. Accordingly, integrating in the unit time causes the cumulative error. If the positioning starts from a beginning point and returns to the same point, the secondary position is not precisely the same initial position and this occurs just because of integration over time. Therefore, the gyroscope's observations are only reliable in the short term, and in the long run there is an error. The complementary filter provides the better of two different accelerometers and gyroscope space. It uses gyroscope observations in short-term because it is very accurate in the short term and is not affected by noise due to the external force input to the platform. And in the long-term observations, this filter only uses accelerometer observations that are much more accurate than the gyroscope's observations and it does not have drift error. As shown in Fig. 1, the integrated gyroscope data in each time epoch would combine with the low-pass data of the accelerometer. The coefficients of 0.98 and 0.02 were used to create the high-pass filter for the gyroscope and the low-pass filter for accelerometers. 2
angle 0.98 (angle gyroData dt ) 0.02 (accData) Fig.1: implementation of the complementary filter [25]
The complementary filter is designed to be repeated in an infinite loop. In each repetition, the roll and pitch angles would be updated by integrating new gyroscope values, then the filter would check that whether the magnitude of acceleration is a reasonable value or not. If the measured acceleration is smaller or larger than the reasonable value, it would be eliminated from computing. 98% of the current the accelerometer observation would be updated with 2% of the new accelerometer observations. By carrying out these actions, it would be ensured that the calculations performed are free of dirt and are very precise [18]. B. Mahony Filter Recently, Mahony, have provided a filter to determine the orientation of a platform, defined in Euclidean 3D space [15]. They have achieved a filter called explicit complement filter or nonlinear complementary filter through changes in the formulation of complementary linear filters and a passive complement [19]. Nonlinear complementary filters include orientation estimators that are inspired by the linear complementary filter. Similar to the linear complementary filter, the non-linear complementary filter is also derived from the combination of vector data of low-frequency orientation (accelerometer and magnetometer) and high-frequency orientation data (gyroscope) [20]. In order to estimate
qE (orientation
E
E
in earth coordinate system), it is necessary to introduce reference vectors a and m . In no-noise mode, the relation between E a and S a is given in equation (1). S
aq q1 E aq q
(1)
In equation (1), the sign , represents the quaternion multiplication and observation taken in the sensor frame, it can be display as , which is displayed as
E
S
aq [0
S
S
ax
mq q1 E mq q
In equation (2), displayed in
S
S
S
ay
a [0 0 g ] in the static mode and g 9.8m.s 2
surface. In the absence of noise and magnetic deviations, the relation between S
aq is the quaternion mode of accelerometers
E
S
az ] . E aq Is the quaternion mode of E a
is the gravitational acceleration on the earth
m and S m is shown as follow.
(2)
mq is the quaternion mode of the magnetometer observations taken in the sensor frame and can be
mq [0
S
mx
S
my
S
mz ] . E mq is the quaternion mode of E m . E m in the absence of magnetic deviation
can be calculated from the world magnetic model (WWM). Equation (3) shows the well-known equation of rigid body
3
motion, which is derived from angular velocity observations of a gyroscope and describes the orientation change in quaternion unit.
q
1 q Sq 2
(3)
S
is the angular velocity resulting from the gyroscope's observations and Sq is the quaternion unit display of S .
S
aq,t qˆt11 E aq,t qˆt 1
(4)
S
mq,t qˆt11 E mq,t qˆt 1
(5)
mes,t S at S aˆt S mt S mˆ t
(6)
ˆb,t ki Smes,t
(7)
S
S
S
ˆ r ,q,t Sq,t 0
k p and ki
S
ˆb,t 0 k p
S
mes,t
(8)
are respectively integral and proportional gains.
1 qˆ qˆt 1 Sˆ r ,q,t 2
(9)
This filter is expressed in various sources called the Explicit Complementary Filter and it is designed for air platforms. The presented algorithm calculates the error by using the multiplication technique of estimated observation of acceleration and magnetic fields, and this error is used to correct the gyroscope bias [21]. C. Madgwick Filter The estimated Euler angles of sensor's frame relative to the Earth's frame are obtained by combining the computational Euler angles from the gyroscope observations and computational Euler angles derived from the combined observations of the accelerometer and the magnetometer measurements. This Integration algorithm is formulated as follows: S ˆest ,t Eq
f t (t ) ( ES qˆest ,t 1 ES q ,t t ) t f
Refer to (10), coefficient
S E
(10)
qˆest ,t and ES qˆest ,t 1 are computational quaternion in the epochs t and t-1 and t is the sampling distance. The
is equal to:
4
t
S E q ,t
t, 1
(11)
On the other hand, the f is equal to:
f J g ,Tb ( ES qˆest ,t 1 , Ebˆ ) f g ,b ( ES qˆest ,t 1 , Saˆ t , Ebˆ , Smˆ t ) t
J gT ( S qˆ E est ,t 1 ) ) T S t J g ( E qˆest ,t 1, Ebˆ ) t
J g ,b (SE qˆest ,t 1, Ebˆ f
(12)
t
f g (SE qˆest ,t 1 , Saˆ ) t ( qˆest ,t 1 , Saˆt , Ebˆ , Smˆ t ) T S t ˆ J g ( E qest ,t 1 , Ebˆ , Smˆ t ) t
S g ,b E
(13)
(14)
In which S aˆt is the accelerometer observation in epoch t, S mˆ is the magnetometer observation in epoch t, and Ebˆ is the t
magnetometer's bias.
Saˆ 0,a ,a
y ,az
(15)
x ,my ,mz
(16)
Ebˆ 0, h 2x h 2y ,0, h z t
(17)
bx h 2x h 2y
(18)
bz hz
(19)
Ehˆ 0, hx , hy , hz ES qˆest ,t 1 Smˆ t ES qˆ*est ,t 1 t
(20)
t
x
Smˆ 0,m t
And also:
J g ( ES qˆest ,t 1 )
2q3 2q2 0
2q4 2q1 2q2 2q1 2q4 2q3 4q2 4q3 0
(21)
J bT ( ES qˆest ,t 1 , Ebˆ ) t
2bz q3 2bz q4 4bx q3 2bz q1 4bx q4 2bz q2
2bz q4 2bz q2 2bx q3 2bx q3 2bz q1 2bx q4 4bz q2 (22) 5 2bx q2 2bz q4 2bx q1 4bz q3 2bx q1 2bz q3 2bx q2
t
fb ( ES qˆest ,t 1 , Ebˆ , Smˆ t ) t
2 bx (0.5 q q 42 ) 2 bz (q2 q4 q1q3 ) mx 2 bx (q2 q3 q1 q4 ) 2 bz (q1 q2 q3 q4 ) m y 2 b (q q q q ) 2 b (0.5 q 2 q 2 ) m z z 2 4 2 3 x 1 3 2 3
(23)
2 (q q q q ) a x 2 4 1 3 (24) f g ( ES qˆest ,t 1 , Saˆt ) 2 (q1 q2 q3 q4 ) a y 2 ( 1 q 2 q 2 ) a z 2 3 2 S Refer to (10), E q ,t is the quaternion derivative or rotational rate of change of the earth frame relative to the sensor frame, which is obtained using the following equation: S E q ,t
1S qˆ Sc,t 2 E est ,t 1
(25)
Refer to (26), Sc ,t is the corrected observations of the gyroscope at time t, which is obtained from the differential gyroscope bias ( Swb,t ) relative to the angular rate observations of a three-axis gyroscope ( S w ). t
Swc,t
Swt
Sw
(26)
b,t
S 0, x , y , z
(27)
Sb,t S ,t t
(28)
t
Refer to (28), the bias of the gyroscope would drift over time with temperature and motion, and this type of error has an accumulation effect in time, while its first value is assumed to be zero. The estimated bias in each epoch is also obtained from the following formula:
Sw ,t 2 ES qˆ*est ,t 1
f f
(29)
Which if we suppose: S ˆest ,t 1 Eq
q1, q2 , q3 , q4
(30)
Then: 6
S * ˆ est ,t 1 Eq
SE qˆest ,t 1 q1, q2 , q3 , q4
Refer to (31),
qˆ
S * E est ,t 1
(31)
is the conjugated quaternion of ES qˆest ,t 1 . is the filter gain, which has the role of displaying the
covariance rate to eliminate measurement errors of the gyroscope with a non-zero mean. On the other hand, refer to (10), also is a filter gain and is called a quaternion derivative of magnitude, which represents all the errors in the mean zero measurement (including sensor noise, similarity of signals to each other in such a way as to be indistinguishable, gradual errors, calibration errors, and non-vertical orientation of sensor axes).To determine and , it is best to use the angle values ω and
respectively. ˆ ˆ 1,ˆ 2 ,ˆ 3 Indicates the static error estimated from the gyroscope's observations for each axis,
and ˆ ˆ 1,ˆ 2 ,ˆ 3 represents the estimated drift rate of the gyroscope bias for each axis. As a result, we have:
1 S 3 ˆest ,t 1 0, ˆ 1 , ˆ 2 , ˆ 3 ˆ Eq 2 4
(32)
1 S 3 ˆest ,t 1 0, ˆ 1 , ˆ 2 , ˆ 3 ˆ 3 Eq 2 4
(33)
By presenting the proposed algorithms for Euler angles determination, the evaluation of its capabilities is discussed in the next section. D. XsensKalman Filter (XKF3i) The Eulerian angles of IMU Sensor (Xsens MTi-G-700)is calculated by Xsens Kalman Filter.XKF3i uses gyroscope, accelerometer and magnetometer observations to calculate the 3D orientation of IMU sensor in which both static and dynamic situations is free from drift error. XKF3i is a proven fusion algorithm that can be found in most Xsens company’s products. The XKF3i fusion algorithm is designed to use the accelerometer (gravity acceleration) and magnetometer (Earth magnetic north) observations to compensate the effect of the drift error resulting from the integration of gyroscope observations (angular velocity from the rate gyros).This kind of drift error compensation is called Attitude and Heading referencing and the system is used to compensate drift error is called Attitude and Heading Reference System (AHRS) [22]. .
3. Results In order to evaluate the accuracy of the presented algorithms and display its capabilities, a field test was conducted on June 11, 2016 in Kish Island. As shown in Fig. 2, a lightweight buoy equipped with an Xsens inertial sensor was tested in this field test and about 25 minutes, test data was collected but just 440 seconds of row data has been calculated by the algorithm as the sample data.
7
Fig. 2. Lightweight buoy used in the test
Fig. 3.Xsens IMU Used in field test
The important point in collecting data is the 8 Hz sampling rate of the inertial sensor. During the Field test, in addition to collecting the Euler angles data, various sensor data including accelerometer, gyroscope and magnetometer were recorded. The observations of these sensors are depicted in Fig. 4 (tri-axis accelerometer observations), Fig. 5 (tri-axis gyroscope observations), and Fig. 6 (tri-axis magnetometer observations).
Fig. 4. Tri-axis accelerometer observations during field testing
8
Fig. 5. Tri-axis gyroscope observations during field testing
Fig. 6. Tri-axis magnetometer observations during field testing
With one view to these figures, it can be inferred that, regardless of noise, the raw observations of these sensors are impregnated with drifts and biases. Using the algorithms and the row data of the sensors, the Euler angles of buoy was determined during the test. Fig. 7, Fig. 8 and Fig. 9 show the calculated Euler angles and the Euler angles obtained from the sensor (XKF3i), respectively.
Fig. 7. Computational(continuous lines) and observational (dotted line) roll angles
9
Fig. 8. Computational(continuous lines) and observational (dotted line) pitch angles
Fig. 9. Computational(continuous lines) and observational (dotted line) yaw angles
In order to show the accuracy of each of the proposed algorithms with reference algorithm (XKF3i) is required to examine the differences betweenXKF3i and each of the algorithms (Complementary, Mahony and Madgwick). Therefore, Fig. 10, Fig. 11 and Fig. 12 show the differences between the proposed algorithms and the XKF3i algorithm to determine the angles of roll, pitch and yaw, respectively.
Fig. 10. Roll angles comparisons between proposed algorithms and XKF3i
10
Fig. 11. Pitch angles comparisons between proposed algorithms and XKF3i
Fig. 12. Yaw angles comparisons between proposed algorithms and XKF3i
Also, the statistical information of these comparisons are presented in TABLE I. Also, the statistical information such as Minimum, Maximum, Standard Deviation and Root Mean Square Deviation that present the accuracy difference between the methods are presented in TABLE I. According to TABLE I, the standard deviation of differences for Madgwick is about 0.57°, 0.37°, 0.50° and also these values for Mahony are 0.56°, 0.39°, 0.50° and for complementary are 0.63°, 0.26°, 2.38° for roll, pitch, and yaw angles, respectively. TABLE I STATISTICAL
INFORMATION
FOR
COMPARISON
BETWEEN
EULERIAN
ANGLES
Madgwick_min Madgwick_max Madgwick_mean Madgwick_std Madgwick_RMSD Mahony_min Mahony_max Mahony_mean Mahony_std Mahony_RMSD Complementary_min Complementary_max Complementary_mean Complementary_std
Roll -0.90 1.23 4.76e-16 0.57 0.57 -1.19 1.20 -1.45e-16 0.56 0.56 -1.12 1.45 -1.74e-16 0.63
11
Pitch -0.82 0.78 1.53e-16 0.37 0.37 -0.89 0.70 8.23e-16 0.39 0.39 -0.70 0.45 7.9e-14 0.26
Yaw -1.15 1.61 4e-16 0.50 0.50 -2.09 1.47 4e-16 0.50 0.50 -3.18 4.88 0.18 2.38
0.63
Complementary_RMSD
0.26
2.38
Standard Deviation Angle (deg)
2.5 2 1.5 1 0.5 0
roll
pitch
yaw
Madg_std
0.57
0.37
0.5
Mah_std
0.56
0.39
0.5
Com_std
0.63
0.26
2.38
Madg_std
Mah_std
Com_std
Fig. 13. Comparison of standard deviation of proposed algorithms
4. Conclusion In this research, Eulerian angles determination algorithms for a buoy was presented using a combination of accelerometer, gyroscope and magnetometer measurements. Then, in order to evaluate the capabilities of this algorithms, a field test was carried out in Kish Island by installing an IMU sensor on a lightweight buoy. According to the statistical information for comparison between Eulerian angles of the proposed algorithms and Eulerian angles of XKF3i, the roll, pitch, and yaw angles differences mean for Madgwick and Mahony algorithms indicate that there is no static bias. On the other hand, the low standard deviations of the Madgwick and Mahony algorithms show the good performance of these algorithms for reducing errors and variable biases with time. By comparing the standard deviation of three algorithms, it is concluded that, Mahony and Madgwick algorithms for determination of roll, pitch and yaw angles have almost equal precision. But if you look a bit more closely at this issue, it would be found that the Mahony algorithm is more accurate than Madgwick algorithm for the roll and yaw angles determination. The results also show that the complementary algorithm is more accurate to determine the pitch angle than the other two algorithms, but is not suitable for determining the yaw angle due to the high standard deviation (high error).Accordingly, it is suggested that appropriate methods for determining particular angles be used in various applications in the field of marine engineering sciences, hydrography and oceanography. APPENDIX Quaternion is a four dimensional complex number used to represent the position of a moving platform or coordinate frame in three dimensional space. An arbitrary condition of frame B relative to frame A can be achieved using the rotation around the axis defined in frame A. As shown in Fig. 8, the axes { Z A ,YA , X A } and { Z B ,YB , X B } Are respectively the frame axis A and B. Refer to (25), { rz , ry , rx } is a unit vector of along the axes in frame A. Quatrain's calculations are often used from the first stage of the normalized quaternion. On this basis, this causes the quaternion units be meter [23]. A B
qˆ q1 , q2 , q3 , q4
[cos , rx sin , ry sin , rz sin ] 2 2 2 2
(25)
12
The conjugate quaternion is shown with the * that is capable of changing the ratio of a frame to another frame. For example, A ˆ is conjugate of AB qˆ and indicates the status of frame A relative to frame B. The conjugate is shown in the equation Bq below: A * ˆ Bq
ABqˆ q1, q2 , q3 , q4
(26)
The quaternion product, denoted by,⊗ can be used to define compound orientations. For example, for two orientations described by BA qˆ and CB qˆ the compounded orientation CA qˆ can be defined by equation: A ˆ CBqˆ BAqˆ Cq
(27)
Fig. 8. Display the status of frame B by a rotation equal to θ around
in
frame A
For two quaternions, a and b, the quaternion product can be determined using the Hamilton rule and defined as (28). A quaternion product is not commutative; that is a b b a .
a b a1 , a2 , a3 , a4 b1 , b2 , b3 , b4 a1b1 a2b2 a3b3 a4b4 a b a2b1 a3b4 a4b3 1 2 a1b3 a2b4 a3b1 a4b2 a1b4 a2b3 a3b2 a4b1
T
(28)
Av
Bv BA qˆ Av BA qˆ*
(29) A B
qˆ
A B
13
R
Bv
A B
R
2q12 1 2q22 2 q2 q3 q1q4 2 q2 q4 q1q3 (30) 2 2 2 q2 q3 q1q4 2q1 1 2q3 2 q3 q4 q1q2 2 2 2 q2 q4 q1q3 2 q3 q4 q1q2 2q1 1 2q4
, ,
a tan 2 2q2 q3 2q1q4 ,2q12 2q22 1
(30)
sin 1 2q2 q4 2q1q3
(31)
tan 1 {2q3q4 2q1q2 }/{2q12 2q42 1}
(32)
References
14
15
Highlights:
Using the commercial algorithm (Xsens Kalman Filter) as the reference for comparing Using magnetometer observations to calculate more precisely in Madwick and Mahony Evaluation of algorithms in the presence or absence of magnetometer observations High sensitivity testing algorithms due to gathering data in non-wavy sea mode
16
S0263-2241(18)31115-1 https://doi.org/10.1016/j.measurement.2018.11.061 MEASUR 6104
To appear in:
Measurement
Received Date: Revised Date: Accepted Date:
23 May 2018 30 October 2018 18 November 2018
Please cite this article as: A. Joybari, H. Amiri, A.A. Ardalan, N.K. Zahraee, Methods Comparison for Attitude Determination of a Lightweight Buoy by Raw Data of IMU, Measurement (2018), doi: https://doi.org/10.1016/ j.measurement.2018.11.061
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Methods Comparison for Attitude Determination of a Lightweight Buoy by Raw Data of IMU Arash Joybari1, Hadi Amiri2, Alireza A. Ardalan3, Niloofar K. Zahraee4 1
Faculty of Engineering and Sustainable Development, University of Gävle, Sweden; [email protected]. School of Engineering Science, College of Engineering, University of Tehran, Iran; [email protected] 3 School of Surveying and Geomatics Engineering, University of Tehran, Iran; [email protected] 3 School of Surveying and Geomatics Engineering, University of Tehran, Iran; [email protected] 2
Correspondence to: Hadi Amiri ([email protected]) Abstract. Today, one of the most important issues is the determination of instantaneous sea level and distinguishing the Tsunami by floating buoy in the ocean. Usually, gyroscopes are used to measure the angular velocity of a buoy. On the other hand, considering the advancement of various technologies in the field of precise accelerometers, make it possible to use these kinds of sensors for navigation purpose. In this research, stable and optimal methods for determining the orientation of a moving buoy is presented using a combination of the gyroscope, accelerometers, and magnetic sensors data. In order to prove the effectiveness of the proposed methods, the raw data were collected from accelerometers, gyroscopes, and magnetometers of (Xsens MTI-G-700) mounted on a Buoy in coastal waters of Kish Island, Iran. Then, by using the proposed methods, the Euler angles of the buoy are determined, while the Euler angles are derived from the Xsens sensor we are considered as a reference. Based on the results, RMSD for Madgwick algorithm are 0.57°, 0.37°, and 0.50° for Mahony algorithm are 0.56°, 0.37°, and 0.50° and finally for Complementary algorithm is 0.63°, 0.26°, and 2.38° which these values are for roll, pitch, and yaw angles respectively. Thus Mahony algorithm for determining roll and yaw Euler angles is more accurate than other algorithms; however, this differences is negligible compared to the Madgwick algorithm. The Complementary algorithm is less accurate than the other two algorithms, especially for determining the yaw angle of the buoy.
1. Introduction Precise measurement of the angles plays a fundamental role in various applications such as aerospace [1], [2], robotics [3], [4], navigation [5], [6], and human movement analysis [7], [8]. Despite the fact that today there are various technologies for measuring the Euler angles, but the use of inertial systems has become more and more progressive since these systems do not depend on the dynamics of movement and environmental conditions [9]. Indeed, a magnetic and inertial measurement unit (MIMU) consists of gyroscope, accelerometer, and magnetometer which, in the case of an orthogonal installation of sensors along the main axis, has the ability to measure rotational and transient Indeed, an IMU (Inertial Measurement Unit) consists of gyroscope, accelerometer, and magnetometer which, in the case of an orthogonal installation of sensors along the main axis, has the ability to measure rotational and transient movements along the three main axes of the body [5]. A gyroscope is able to measure the sensor’s attitude during the time, if its initial state is known [10], [11]. However, because of 1
the accumulation errors in gyroscope measurements, the calculated Euler angles can have significant errors. Accordingly, accelerometers and magnetometers, which measure gravity and magnetic fields, can be integrated with gyroscopes and thus produce accurate Euler angles. It is worth noting that accelerometers are also exposed to relatively high noise levels, such that sudden motion of the sensor leads to a reduction in the accuracy of the accelerometer's gravity observations. In general, the use of an appropriate integration algorithm for calculating and estimating optimal Euler angles is essential by combining gyroscope, accelerometer, and magnetometer measurements [12]. Typically, the Kalman filter has been used as a wellknown tool in most of the commonly used research in Euler angles determination using raw observations of gyroscope, accelerometer, and magnetometer [13]. The basis of the Kalman process is based on the repetition of linear or non-linear regressions in different epochs. The widespread use of the Kalman filter is evidence of its high precision and performance, although it has some disadvantages. Alternatively, we can use the fuzzy process or fixed filters to determine the Euler angles in order to improve the accuracy of accelerometer observations at low angular velocities and gyroscope observations at high angular velocities [14]. On the other hand, [15] proposed another method based on least squares regression, so that the filter integrates the observations of all inertial sensors in determining the Euler angles. In this research, three methods for integrating the inertial sensor observations are presented in determining the Euler angles while in our previous research just two methods for integrating the inertial sensor observations has been evaluated [16]. The fundamental of this filter is based on the use of quaternion for moving platform angles determination in a three-dimensional space [17]. The Comparison of proposed filters and commercial filters (Xsens Kalman filter) has been investigated. In the next section, the details of the proposed method are discussed. After that, the results of applying the proposed method are presented.
2. Materials and Methods A. Complementary Filter The Kalman filter is very efficient and accurate, but it also has some disadvantages as complexity. The complementary filter is much easier to understand and implementation than the Kalman filter. Most IMUs have 6 degrees of freedom, which means that there are 3 accelerometers and 3 gyroscopes. The IMU has the ability to measure the position and orientation of the platform to which it is attached because an object in space has 6 degrees of freedom. So, if 6 degrees of freedom for an object in space is measured, all problems related to determining of position and orientation of object would be solved. The raw data of the sensors in the IMU unit is not enough to achieve this goal. Observations of accelerometer and gyroscope sensors are used to reach the position of a platform in space. The gyroscope's observations are angular velocities; therefore, by integrating the angular velocity per time unit, the Eulerian angles are obtained. The accelerometer observations calculate the Eulerian angles, though. Accelerometers measure the position and state of the gravity vector (g-force) and calculate the platform orientation using the trigonometric relationship. A major problem with using these two types of observation to reach the platform orientation is that it is not possible to determine the orientation of platform without filtering on them. The problem with accelerometer measurements is that the accelerometers, in addition to gravitational acceleration, measure all the forces on the moving platform. The accelerometer, in addition to measuring gravitational acceleration, also measures noises, so accelerometer observations are reliable only in long-term measurements and require a low pass filter for removing noises. But the gyroscope calculates the Eulerian angles by integrating the angular rate in the unit time. Accordingly, integrating in the unit time causes the cumulative error. If the positioning starts from a beginning point and returns to the same point, the secondary position is not precisely the same initial position and this occurs just because of integration over time. Therefore, the gyroscope's observations are only reliable in the short term, and in the long run there is an error. The complementary filter provides the better of two different accelerometers and gyroscope space. It uses gyroscope observations in short-term because it is very accurate in the short term and is not affected by noise due to the external force input to the platform. And in the long-term observations, this filter only uses accelerometer observations that are much more accurate than the gyroscope's observations and it does not have drift error. As shown in Fig. 1, the integrated gyroscope data in each time epoch would combine with the low-pass data of the accelerometer. The coefficients of 0.98 and 0.02 were used to create the high-pass filter for the gyroscope and the low-pass filter for accelerometers. 2
angle 0.98 (angle gyroData dt ) 0.02 (accData) Fig.1: implementation of the complementary filter [25]
The complementary filter is designed to be repeated in an infinite loop. In each repetition, the roll and pitch angles would be updated by integrating new gyroscope values, then the filter would check that whether the magnitude of acceleration is a reasonable value or not. If the measured acceleration is smaller or larger than the reasonable value, it would be eliminated from computing. 98% of the current the accelerometer observation would be updated with 2% of the new accelerometer observations. By carrying out these actions, it would be ensured that the calculations performed are free of dirt and are very precise [18]. B. Mahony Filter Recently, Mahony, have provided a filter to determine the orientation of a platform, defined in Euclidean 3D space [15]. They have achieved a filter called explicit complement filter or nonlinear complementary filter through changes in the formulation of complementary linear filters and a passive complement [19]. Nonlinear complementary filters include orientation estimators that are inspired by the linear complementary filter. Similar to the linear complementary filter, the non-linear complementary filter is also derived from the combination of vector data of low-frequency orientation (accelerometer and magnetometer) and high-frequency orientation data (gyroscope) [20]. In order to estimate
qE (orientation
E
E
in earth coordinate system), it is necessary to introduce reference vectors a and m . In no-noise mode, the relation between E a and S a is given in equation (1). S
aq q1 E aq q
(1)
In equation (1), the sign , represents the quaternion multiplication and observation taken in the sensor frame, it can be display as , which is displayed as
E
S
aq [0
S
S
ax
mq q1 E mq q
In equation (2), displayed in
S
S
S
ay
a [0 0 g ] in the static mode and g 9.8m.s 2
surface. In the absence of noise and magnetic deviations, the relation between S
aq is the quaternion mode of accelerometers
E
S
az ] . E aq Is the quaternion mode of E a
is the gravitational acceleration on the earth
m and S m is shown as follow.
(2)
mq is the quaternion mode of the magnetometer observations taken in the sensor frame and can be
mq [0
S
mx
S
my
S
mz ] . E mq is the quaternion mode of E m . E m in the absence of magnetic deviation
can be calculated from the world magnetic model (WWM). Equation (3) shows the well-known equation of rigid body
3
motion, which is derived from angular velocity observations of a gyroscope and describes the orientation change in quaternion unit.
q
1 q Sq 2
(3)
S
is the angular velocity resulting from the gyroscope's observations and Sq is the quaternion unit display of S .
S
aq,t qˆt11 E aq,t qˆt 1
(4)
S
mq,t qˆt11 E mq,t qˆt 1
(5)
mes,t S at S aˆt S mt S mˆ t
(6)
ˆb,t ki Smes,t
(7)
S
S
S
ˆ r ,q,t Sq,t 0
k p and ki
S
ˆb,t 0 k p
S
mes,t
(8)
are respectively integral and proportional gains.
1 qˆ qˆt 1 Sˆ r ,q,t 2
(9)
This filter is expressed in various sources called the Explicit Complementary Filter and it is designed for air platforms. The presented algorithm calculates the error by using the multiplication technique of estimated observation of acceleration and magnetic fields, and this error is used to correct the gyroscope bias [21]. C. Madgwick Filter The estimated Euler angles of sensor's frame relative to the Earth's frame are obtained by combining the computational Euler angles from the gyroscope observations and computational Euler angles derived from the combined observations of the accelerometer and the magnetometer measurements. This Integration algorithm is formulated as follows: S ˆest ,t Eq
f t (t ) ( ES qˆest ,t 1 ES q ,t t ) t f
Refer to (10), coefficient
S E
(10)
qˆest ,t and ES qˆest ,t 1 are computational quaternion in the epochs t and t-1 and t is the sampling distance. The
is equal to:
4
t
S E q ,t
t, 1
(11)
On the other hand, the f is equal to:
f J g ,Tb ( ES qˆest ,t 1 , Ebˆ ) f g ,b ( ES qˆest ,t 1 , Saˆ t , Ebˆ , Smˆ t ) t
J gT ( S qˆ E est ,t 1 ) ) T S t J g ( E qˆest ,t 1, Ebˆ ) t
J g ,b (SE qˆest ,t 1, Ebˆ f
(12)
t
f g (SE qˆest ,t 1 , Saˆ ) t ( qˆest ,t 1 , Saˆt , Ebˆ , Smˆ t ) T S t ˆ J g ( E qest ,t 1 , Ebˆ , Smˆ t ) t
S g ,b E
(13)
(14)
In which S aˆt is the accelerometer observation in epoch t, S mˆ is the magnetometer observation in epoch t, and Ebˆ is the t
magnetometer's bias.
Saˆ 0,a ,a
y ,az
(15)
x ,my ,mz
(16)
Ebˆ 0, h 2x h 2y ,0, h z t
(17)
bx h 2x h 2y
(18)
bz hz
(19)
Ehˆ 0, hx , hy , hz ES qˆest ,t 1 Smˆ t ES qˆ*est ,t 1 t
(20)
t
x
Smˆ 0,m t
And also:
J g ( ES qˆest ,t 1 )
2q3 2q2 0
2q4 2q1 2q2 2q1 2q4 2q3 4q2 4q3 0
(21)
J bT ( ES qˆest ,t 1 , Ebˆ ) t
2bz q3 2bz q4 4bx q3 2bz q1 4bx q4 2bz q2
2bz q4 2bz q2 2bx q3 2bx q3 2bz q1 2bx q4 4bz q2 (22) 5 2bx q2 2bz q4 2bx q1 4bz q3 2bx q1 2bz q3 2bx q2
t
fb ( ES qˆest ,t 1 , Ebˆ , Smˆ t ) t
2 bx (0.5 q q 42 ) 2 bz (q2 q4 q1q3 ) mx 2 bx (q2 q3 q1 q4 ) 2 bz (q1 q2 q3 q4 ) m y 2 b (q q q q ) 2 b (0.5 q 2 q 2 ) m z z 2 4 2 3 x 1 3 2 3
(23)
2 (q q q q ) a x 2 4 1 3 (24) f g ( ES qˆest ,t 1 , Saˆt ) 2 (q1 q2 q3 q4 ) a y 2 ( 1 q 2 q 2 ) a z 2 3 2 S Refer to (10), E q ,t is the quaternion derivative or rotational rate of change of the earth frame relative to the sensor frame, which is obtained using the following equation: S E q ,t
1S qˆ Sc,t 2 E est ,t 1
(25)
Refer to (26), Sc ,t is the corrected observations of the gyroscope at time t, which is obtained from the differential gyroscope bias ( Swb,t ) relative to the angular rate observations of a three-axis gyroscope ( S w ). t
Swc,t
Swt
Sw
(26)
b,t
S 0, x , y , z
(27)
Sb,t S ,t t
(28)
t
Refer to (28), the bias of the gyroscope would drift over time with temperature and motion, and this type of error has an accumulation effect in time, while its first value is assumed to be zero. The estimated bias in each epoch is also obtained from the following formula:
Sw ,t 2 ES qˆ*est ,t 1
f f
(29)
Which if we suppose: S ˆest ,t 1 Eq
q1, q2 , q3 , q4
(30)
Then: 6
S * ˆ est ,t 1 Eq
SE qˆest ,t 1 q1, q2 , q3 , q4
Refer to (31),
qˆ
S * E est ,t 1
(31)
is the conjugated quaternion of ES qˆest ,t 1 . is the filter gain, which has the role of displaying the
covariance rate to eliminate measurement errors of the gyroscope with a non-zero mean. On the other hand, refer to (10), also is a filter gain and is called a quaternion derivative of magnitude, which represents all the errors in the mean zero measurement (including sensor noise, similarity of signals to each other in such a way as to be indistinguishable, gradual errors, calibration errors, and non-vertical orientation of sensor axes).To determine and , it is best to use the angle values ω and
respectively. ˆ ˆ 1,ˆ 2 ,ˆ 3 Indicates the static error estimated from the gyroscope's observations for each axis,
and ˆ ˆ 1,ˆ 2 ,ˆ 3 represents the estimated drift rate of the gyroscope bias for each axis. As a result, we have:
1 S 3 ˆest ,t 1 0, ˆ 1 , ˆ 2 , ˆ 3 ˆ Eq 2 4
(32)
1 S 3 ˆest ,t 1 0, ˆ 1 , ˆ 2 , ˆ 3 ˆ 3 Eq 2 4
(33)
By presenting the proposed algorithms for Euler angles determination, the evaluation of its capabilities is discussed in the next section. D. XsensKalman Filter (XKF3i) The Eulerian angles of IMU Sensor (Xsens MTi-G-700)is calculated by Xsens Kalman Filter.XKF3i uses gyroscope, accelerometer and magnetometer observations to calculate the 3D orientation of IMU sensor in which both static and dynamic situations is free from drift error. XKF3i is a proven fusion algorithm that can be found in most Xsens company’s products. The XKF3i fusion algorithm is designed to use the accelerometer (gravity acceleration) and magnetometer (Earth magnetic north) observations to compensate the effect of the drift error resulting from the integration of gyroscope observations (angular velocity from the rate gyros).This kind of drift error compensation is called Attitude and Heading referencing and the system is used to compensate drift error is called Attitude and Heading Reference System (AHRS) [22]. .
3. Results In order to evaluate the accuracy of the presented algorithms and display its capabilities, a field test was conducted on June 11, 2016 in Kish Island. As shown in Fig. 2, a lightweight buoy equipped with an Xsens inertial sensor was tested in this field test and about 25 minutes, test data was collected but just 440 seconds of row data has been calculated by the algorithm as the sample data.
7
Fig. 2. Lightweight buoy used in the test
Fig. 3.Xsens IMU Used in field test
The important point in collecting data is the 8 Hz sampling rate of the inertial sensor. During the Field test, in addition to collecting the Euler angles data, various sensor data including accelerometer, gyroscope and magnetometer were recorded. The observations of these sensors are depicted in Fig. 4 (tri-axis accelerometer observations), Fig. 5 (tri-axis gyroscope observations), and Fig. 6 (tri-axis magnetometer observations).
Fig. 4. Tri-axis accelerometer observations during field testing
8
Fig. 5. Tri-axis gyroscope observations during field testing
Fig. 6. Tri-axis magnetometer observations during field testing
With one view to these figures, it can be inferred that, regardless of noise, the raw observations of these sensors are impregnated with drifts and biases. Using the algorithms and the row data of the sensors, the Euler angles of buoy was determined during the test. Fig. 7, Fig. 8 and Fig. 9 show the calculated Euler angles and the Euler angles obtained from the sensor (XKF3i), respectively.
Fig. 7. Computational(continuous lines) and observational (dotted line) roll angles
9
Fig. 8. Computational(continuous lines) and observational (dotted line) pitch angles
Fig. 9. Computational(continuous lines) and observational (dotted line) yaw angles
In order to show the accuracy of each of the proposed algorithms with reference algorithm (XKF3i) is required to examine the differences betweenXKF3i and each of the algorithms (Complementary, Mahony and Madgwick). Therefore, Fig. 10, Fig. 11 and Fig. 12 show the differences between the proposed algorithms and the XKF3i algorithm to determine the angles of roll, pitch and yaw, respectively.
Fig. 10. Roll angles comparisons between proposed algorithms and XKF3i
10
Fig. 11. Pitch angles comparisons between proposed algorithms and XKF3i
Fig. 12. Yaw angles comparisons between proposed algorithms and XKF3i
Also, the statistical information of these comparisons are presented in TABLE I. Also, the statistical information such as Minimum, Maximum, Standard Deviation and Root Mean Square Deviation that present the accuracy difference between the methods are presented in TABLE I. According to TABLE I, the standard deviation of differences for Madgwick is about 0.57°, 0.37°, 0.50° and also these values for Mahony are 0.56°, 0.39°, 0.50° and for complementary are 0.63°, 0.26°, 2.38° for roll, pitch, and yaw angles, respectively. TABLE I STATISTICAL
INFORMATION
FOR
COMPARISON
BETWEEN
EULERIAN
ANGLES
Madgwick_min Madgwick_max Madgwick_mean Madgwick_std Madgwick_RMSD Mahony_min Mahony_max Mahony_mean Mahony_std Mahony_RMSD Complementary_min Complementary_max Complementary_mean Complementary_std
Roll -0.90 1.23 4.76e-16 0.57 0.57 -1.19 1.20 -1.45e-16 0.56 0.56 -1.12 1.45 -1.74e-16 0.63
11
Pitch -0.82 0.78 1.53e-16 0.37 0.37 -0.89 0.70 8.23e-16 0.39 0.39 -0.70 0.45 7.9e-14 0.26
Yaw -1.15 1.61 4e-16 0.50 0.50 -2.09 1.47 4e-16 0.50 0.50 -3.18 4.88 0.18 2.38
0.63
Complementary_RMSD
0.26
2.38
Standard Deviation Angle (deg)
2.5 2 1.5 1 0.5 0
roll
pitch
yaw
Madg_std
0.57
0.37
0.5
Mah_std
0.56
0.39
0.5
Com_std
0.63
0.26
2.38
Madg_std
Mah_std
Com_std
Fig. 13. Comparison of standard deviation of proposed algorithms
4. Conclusion In this research, Eulerian angles determination algorithms for a buoy was presented using a combination of accelerometer, gyroscope and magnetometer measurements. Then, in order to evaluate the capabilities of this algorithms, a field test was carried out in Kish Island by installing an IMU sensor on a lightweight buoy. According to the statistical information for comparison between Eulerian angles of the proposed algorithms and Eulerian angles of XKF3i, the roll, pitch, and yaw angles differences mean for Madgwick and Mahony algorithms indicate that there is no static bias. On the other hand, the low standard deviations of the Madgwick and Mahony algorithms show the good performance of these algorithms for reducing errors and variable biases with time. By comparing the standard deviation of three algorithms, it is concluded that, Mahony and Madgwick algorithms for determination of roll, pitch and yaw angles have almost equal precision. But if you look a bit more closely at this issue, it would be found that the Mahony algorithm is more accurate than Madgwick algorithm for the roll and yaw angles determination. The results also show that the complementary algorithm is more accurate to determine the pitch angle than the other two algorithms, but is not suitable for determining the yaw angle due to the high standard deviation (high error).Accordingly, it is suggested that appropriate methods for determining particular angles be used in various applications in the field of marine engineering sciences, hydrography and oceanography. APPENDIX Quaternion is a four dimensional complex number used to represent the position of a moving platform or coordinate frame in three dimensional space. An arbitrary condition of frame B relative to frame A can be achieved using the rotation around the axis defined in frame A. As shown in Fig. 8, the axes { Z A ,YA , X A } and { Z B ,YB , X B } Are respectively the frame axis A and B. Refer to (25), { rz , ry , rx } is a unit vector of along the axes in frame A. Quatrain's calculations are often used from the first stage of the normalized quaternion. On this basis, this causes the quaternion units be meter [23]. A B
qˆ q1 , q2 , q3 , q4
[cos , rx sin , ry sin , rz sin ] 2 2 2 2
(25)
12
The conjugate quaternion is shown with the * that is capable of changing the ratio of a frame to another frame. For example, A ˆ is conjugate of AB qˆ and indicates the status of frame A relative to frame B. The conjugate is shown in the equation Bq below: A * ˆ Bq
ABqˆ q1, q2 , q3 , q4
(26)
The quaternion product, denoted by,⊗ can be used to define compound orientations. For example, for two orientations described by BA qˆ and CB qˆ the compounded orientation CA qˆ can be defined by equation: A ˆ CBqˆ BAqˆ Cq
(27)
Fig. 8. Display the status of frame B by a rotation equal to θ around
in
frame A
For two quaternions, a and b, the quaternion product can be determined using the Hamilton rule and defined as (28). A quaternion product is not commutative; that is a b b a .
a b a1 , a2 , a3 , a4 b1 , b2 , b3 , b4 a1b1 a2b2 a3b3 a4b4 a b a2b1 a3b4 a4b3 1 2 a1b3 a2b4 a3b1 a4b2 a1b4 a2b3 a3b2 a4b1
T
(28)
Av
Bv BA qˆ Av BA qˆ*
(29) A B
qˆ
A B
13
R
Bv
A B
R
2q12 1 2q22 2 q2 q3 q1q4 2 q2 q4 q1q3 (30) 2 2 2 q2 q3 q1q4 2q1 1 2q3 2 q3 q4 q1q2 2 2 2 q2 q4 q1q3 2 q3 q4 q1q2 2q1 1 2q4
, ,
a tan 2 2q2 q3 2q1q4 ,2q12 2q22 1
(30)
sin 1 2q2 q4 2q1q3
(31)
tan 1 {2q3q4 2q1q2 }/{2q12 2q42 1}
(32)
References
14
15
Highlights:
Using the commercial algorithm (Xsens Kalman Filter) as the reference for comparing Using magnetometer observations to calculate more precisely in Madwick and Mahony Evaluation of algorithms in the presence or absence of magnetometer observations High sensitivity testing algorithms due to gathering data in non-wavy sea mode
16