A cross product calibration method for micro-electro mechanical system gyroscope in unmanned aerial vehicle attitude determination system

This paper presents a novel calibration method for micro-electro mechanical system gyroscope in attitude measurement system of small rotor unmanned aerial vehicles. This method is based on an observation vector and its cross product, which is especially valuable for the in-field calibration without the aid of external equipment. By analysing the error model of the tri-axial gyroscope, the principle of calibration is proposed. Compared with other algorithms, numerical simulations are performed to evaluate the effectiveness of integral form of the cross product calibration method. Experiment on the hex-rotor unmanned aerial vehicle platform shows that the proposed method has great advantages in low-cost integrated navigation system.


Introduction
Attitude determination of unmanned aerial vehicles (UAVs) is essentially an integrated navigation problem. 1 The theory and technique of the calibration for conventional inertial devices (especially in integrated navigation systems (INS)) are well established. [2][3][4][5] These calibration methods all rely on high-precision equipment, which can provide accurate attitude reference and angular velocity. [6][7][8][9][10] For an INS equipped with low-end micro-electro mechanical system (MEMS) sensors, using specific and precision equipment to calibrate tri-axial gyroscopes will raise the cost and is usually impractical. Therefore, the above-mentioned calibration methods are probably unsuitable, and using attitude information for calibration is the best strategy.
Liu and Fang proposed a six-position rotation calibration method for field calibration of a fibre optic gyroscope (FOG) inertial measurement unit (IMU). 11 The reported experimental results show that the calibration accuracy is at the same level as the conventional method using a turntable. However, it is noteworthy that the three axes of the IMU are required to point to the sky and the earth in turn and rotate 360 degrees in both directions. If there are errors in the actual direction or rotation angle, the calibration results will be directly affected. Solutions to this problem are not described in detail in their research reports.
Wan and Huang proposed a gyroscope field calibration algorithm based on heading angle, pitch angle and roll angle information. The pitch angle and the roll angle were measured by an accelerometer, and the heading angle was obtained by a magnetic compass or the Global Positioning System. 12 W. T. Fong et al. proposed a method for IMU calibration without external equipment. This method directly refers to gravitational acceleration and realizes the calibration of tri-axial gyroscope. 13 Cheuk et al. 14 used the W. T. Fong's method for IMU self-calibration. E. Dorveaux et al. proposed another self-calibration method of a tri-axial gyroscope with reference to gravitational acceleration. [15][16][17][18][19][20] According to the research of X. Li and Z. Li, 21 a novel calibration method that employs the invariance of the dot product of two constant vectors is presented for tri-axial field sensors in strap-down navigation systems, and this method has been used to calibrate threeaxis magnetometer. But we cannot find an auxiliary vector that has constant dot product with the angular velocity vector, so this method is not suitable for triaxial gyroscope.
Based on the above analysis, a vector cross product calibration method for a tri-axial gyroscope is proposed, which is expressed as the cross product of the reference vector and the angular velocity. This method is derived from the derivative of the reference vector with respect to time, which can be used in INS.

Methodology
Attitude algorithms for INS rely on three 3D vectors, namely the geomagnetic vector, gravity vector and angular velocity. In ideal cases, the sensor output should be proportional to the measurand, that is, v = ku, in which k is the sensitivity (scaling factor), cate the vector being measured and the sensor output, respectively. However, the sensor usually has different error sources that come from manufacture and fabrication, such as bias, scale factor, misalignment and noise. 21,22 A generalized error model of the tri-axial vector sensor can be described by equation (1 In equation (1), the vector e = ( e 1 e 2 e 3 ) T represents sensor noise and higher order error terms, and it can be further divided into the sensor noise e 0 and the higher order terms (i.e. the non-linear part) H.O.T. The calibration of low-end three-axis gyroscope is usually aimed at constant error and linear error, that is, K and b in equation (1), while the terms of e 0 and higher order error are neglected. Thus, equation (1) will become a linear error model as v = Ku + b, in which all the constant errors can be included in b = ( b 1 b 2 b 3 ) T , and all the errors that are proportional to vector u can be described by the matrix K = (k ij ) 333 .
For a low-end tri-axial gyroscope, if the angular velocity in the aircraft coordinate frame is v, while the output of gyroscope isṽ, the error model can be simplified as equation ( For a time-invariant vector u in the geographic coordinate frame, its time derivative can be calculated by its cross product with v, as shown in equation (3) Equation (3) can also be written in the form of matrix multiplication, as shown in equation (4) In equation (4), the matrix ½v3 on the right is called 'cross product matrix', which corresponds to v. For any non-zero vector u, u½v3 = u3v if u is written as a row vector, and ½v3u = v3u if u is written as a column vector. If K and b are determined, we denote L = K À1 and d = K À1 b, then substitute the error model equation (2) into equation (4), and it leads to equation (5) According to equation (5), if u andṽ are both available, Land d can be determined, and the calibration of gyroscope can be implemented. Equation (5) is the basis of the vector-aided calibration method.
In MEMS-based INS, the vector u can be provided by calibrated tri-axial magnetometer or accelerometer, that is, the geomagnetic or gravity vector can serve as u.
In the process of data acquisition, the INS needs to be rotated to provide an excitation signal, that is, a nonzero angular velocity.
Equation (5) is in a differential form. Naturally we can transform equation (5) into an integral form. Assuming that the starting and ending time of the rotation process are t 0 and T, respectively, the change of vector u (denoted as Du) can be calculated by equation (6 To get a closed-form solution of the elements in L and d, the symbol I(u iṽj ) is introduced to represent the integral I(u iṽj )[ R T t 0 u iṽj dt, and equation (6) can also be rewritten in a component-wise form, which is shown in equation (7) Du 1 Equation (7) is a linear expression of all the elements in L and d, and thus it is a linear regression problem.
Equations (5) and (6) are both derived from the cross product in equation (3), and the only difference is that equation (5) is in differential form, while equation (6) is in integral form. We will focus on the integral form in the following discussion.
It can be seen that the above-mentioned cross product calibration method shares the same principle with the calibration method described in previous studies, [13][14][15][16][17][18][19][20][21] but the proposed method greatly simplifies the calculation in the calibration process. Meanwhile, the proposed calibration method only needs a time-invariant vector u (gravity vector g or geomagnetic vector h) without using turntable or other precision equipment. Moreover, when using the magnetometer or accelerometer to provide the reference vector u in calibration, the misalignment between gyroscope and magnetometer or accelerometer can be eliminated simultaneously.

Numerical simulation and experiment
In the following simulations, we use the North-East-Down (NED) coordinate frame and use the gravity vector g 0 = (0 0 9.8) T m s 22 as the reference vector u.
We use the trapezoidal method (described in Sa¨rkkaë t al. 23 ) to calculate the integration in equation (7) numerically. The sampling rate is 100 Hz, and the corresponding time interval is 0.01 s.

Rotation arrangement in calibration
To implement the gyro calibration, we first collect raw data during several rotations of the UAV. First, we point the x, y and z axis of the IMU vertically upwards and downwards, respectively, and thus there are six different cases. Then, we rotate the IMU both clockwise and counter-clockwise around x, y and z axis, respectively, so that a total of 36 different rotations can be obtained. Furthermore, in each rotation, the angular velocity increases linearly from 0 to 90°/s within 1 s and then decreases linearly to zero again within 1 s. So the rotation angle is exactly 90.
According to the above arrangement, there will be 36 rotations corresponding to different combinations of the angular velocity v and the reference vector u. However, among these 36 rotations, there are 12 rotations in which the angular velocity is parallel or antiparallel to the reference vector, and thus the reference vector u is unchanged during the rotation. In such cases, both sides of equations (5) and (6) become zero. In other words, these 12 rotations do not provide effective data for gyroscope calibration. In fact, this kind of rotation is called 'useless rotation' in Li et al. 24 In contrast, the reference vector and the angular velocity are orthogonal to each other in the other 24 rotations, and such rotations can be referred to as 'effective rotations', as shown in Figure 1.
As stated in section 'Introduction,' the principle of the proposed cross product calibration method is quite similar to that of the Fong's method. 13 In fact, if the reference vector u on the right side of equation (7) is calculated by differential equation method, equation (7) can be regarded as a first-order approximation of Fong's method. However, in Fong's method, the reference vector u is calculated by INS updating algorithm (such as the fourth-order Runge-Kutta algorithm in the following simulation), rather than provided by accelerometer. As a result, the elements of L and d will be   (7) and that will inevitably lead to deviation of the calibration results. What is more, it can be expected that such deviation will increase with larger gyroscope error.
To compare Fong's calibration method with the proposed algorithm, we carry out the following numerical simulations. For Fong's method, we use the fourthorder Runge-Kutta (R-K) algorithm to calculate the reference vector (i.e. the gravity vector), as shown in equation (8) Equation (8) describes the renewal process of gravity vector from the moment k to k + 1. In equation (8), h is the time interval of the numerical algorithm. It should be noted thatṽ k + 0:5 in equation (8) represents the measured angular velocity at h/2 after the moment k. In other words, the sampling rate of the gyroscope must be twice of that of the reference vector. Therefore, for Fong's method, the sampling rate of gyroscope is set to 200 Hz.
We assume that the gyro has the following errors We first compare the proposed cross product calibration method with Fong's method with different rotation arrangements, that is, the above-mentioned 36 and 24 rotations. The simulation results are presented in Table 1.
According to Table 1, the performance of the cross product calibration method is good with either 36 or 24 rotations, but Fong's method shows unsatisfactory results, especially with 36 rotations. That means Fong's method is more significantly affected by the useless rotations.
Next, we evaluate the impacts of the rotation angle, which is also an important issue. It can be expected that the cumulative error and its influence on Fong's calibration method will increase with larger rotation angle. We set the rotation angle to 180°, 270°and 360°, respectively, and the simulation results are presented in Table 2.
According to Table 2, the result of Fong's method is still poorer than that of the proposed method when the rotation angle is 180°or 270°. However, it is particularly noteworthy that both methods cannot calculate the matrix L correctly when the rotation angle is 360°. This is mainly because the reference vector has no change after a 360°rotation and thus it can be viewed as another kind of 'useless rotation'.
According to the above simulation results, the proposed cross product calibration method can give accurate results in most cases, but the rotation range should be less than 360°. Therefore, when using the proposed cross product calibration method, a sufficient and feasible rotation arrangement is as follows: 1. Point the three axes of the IMU vertically upwards and downwards, respectively (six different cases in total). 2. In each of the above cases, rotate the IMU around the two horizontal axes both clockwise and counter-clockwise (see Figure 1, 24 rotations in total). 3. The angle of each rotation is 90°. We will use this arrangement in the following numerical simulations.

Impacts of sensor noise and error magnitude
To make the test more comprehensive, we randomly select 1000 error samples for the numerical simulation.
The distortion matrix is generated as L = I 333 + D, where I 333 is an identity matrix, and each element of D is zero-mean and normally distributed with the standard deviation s D = 0:01. Besides, the biases are also assumed to be normally distributed with the standard deviation s b = 3 8 =s.
To evaluate the calibration results of the matrix L, we use an evaluation factor based on the Frobenius  matrix norm. If the presumed distortion matrix is L c , and the simulation result is L p , the evaluation factor will be defined as In equation (9), the symbol|Á| F indicates the Frobenius norm of a matrix. For any 333 matrix A, |I2A| F reflects the difference between A and Identity matrix I. Obviously, if L p is close to L c , Ev(L) will be close to zero. On the contrary, if Ev(L) is close to (or even larger than) 1, it means that the calibration method cannot eliminate the error of gyroscope.
On the contrary, we evaluate the calibration result of the bias by the relative error Ev(b), as shown in equation (10) In equation (10), b p and b c denote the presumed bias and the simulation result, respectively.
First, we add a Gaussian white noise item to the measurement of gyroscope, with the standard deviation s v varying from 0.2 to 2.0°/s. The simulation results are shown in Figure 2. As can be seen from Figure 2, when the noise of gyroscope enlarges, the evaluation factors Ev(L) and Ev(b) of the two calibration methods both increase, but the results of the proposed method are much better than those of Fong's method.
Second, we evaluate the calibration results when the reference vector has noise. We add a Gaussian white noise item to the gravity vector, and the standard deviation varies from 1 to 10 mg. The calibration effect is shown in Figure 3. Once again, the proposed method shows better performance than that of Fong's method. Nevertheless, we can also see that the noise of reference vector has greater impacts to the proposed method, since this noise item will affect both sides of equation (7).
Moreover, we test the effectiveness of both calibration methods when dealing with different magnitudes of sensor errors. In this test, s D varies from 0.02 to 0.1, and s b changes from 2.5 to 10°/s, while the noise level s v is set to 0.1°/s. The simulation results are presented in Figures 4 and 5. It is clear that the proposed method always shows better performance compared to Fong's method, even with large s D and s b .

Experiment on UAV
We apply all the above-mentioned methods to a hexrotor UAV, which consists of a tri-axial magnetometer   (HMC5883), a tri-axial accelerometer (ADXL345) and a MEMS gyroscope (L3G4200D), as shown in Figure 6.
In the experiment, the UAV is rotated in hand, and different calibration methods of raw data are collected. Data are collected in 24 rotations following the recommended procedure, that is, approximately 90°rotations around the horizontal axes. The accelerometer ADXL345 is calibrated before the experiment, in order to serve as the reference for the calibration of the gyroscope, and then L and b calculated by different calibration methods are compensated.
Finally, we installed the UAV on the test platform for static live test and observe hovering test, as shown in Figure 7. Figure 8 shows the comparison of gyroscope yaw angular deviations before and after calibration in the static live test. In Figure 8, the blue line represents the yaw angle calculated according to uncalibrated gyroscope data, which shows the drift up to 25.6°in 1 min. The red line represents the yaw angle according to gyroscope data calibrated by Fong's method, and the drift error can keep within 8.0°. The green line represents the error of yaw angle after calibration with cross product algorithm, and the effect is better than that of Fong's method. It can be seen from Figure 8 that the attitude angle would accumulate drift after the solution, because there still have been some errors after the calibration of the UAV gyroscope. The corrected sampling value can be combined with the data fusion algorithm, 25,26 for example, Kalman filter and complementary filter, to suppress the accumulated errors and ensure the attitude stability for long-time flight.
Next, the UAV remote control throttle remains 50% for hovering test. Figure 9 shows the comparison of yaw angular deviations after calibration with two methods. As we can see in Figure 9, the red line represents the yaw angle obtained by Fong's method, and the drift error of yaw angle can remain within 1.5°. The green   line represents the error of yaw angle after calibration with the cross product algorithm, and the effect is obviously better than that of Fong's method.

Conclusion
We introduce a gravity-aided calibration method for low-end MEMS gyros in this paper. The proposed method essentially unifies the calibration method based on gravity vector, which makes its expression more concise and greatly simplifies the calculation in the calibration process. This method requires no external equipment and is particularly suitable for the in-field calibration of a MEMS gyroscope in a UAV. Although its precision may be limited in practice, this method is efficient to provide a good initial value for the online parameter estimation. The proposed method is based on the time derivative of the gravity vector, and the integral form cross product is used. Numerical simulations are carried out to test the performance of the cross product calibration method and to imitate Fong's method. In the experiment on a UAV platform, the results show that the proposed method can offer higher accuracy in most cases.
The basis of the proposed method, that is, the time derivative of the gravity vector, can be replaced by the time derivative of any constant vector in the aircraft coordinate frame. Therefore, such a vector can be utilized instead of the gravity, and the proposed method will then be evolved to a vector-aided scheme.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.