A Novel MEMS Gyroscope In-Self Calibration Approach

This paper presents a novel approach for hand-held low-cost MEMS (micro-electro-mechanical system) gyroscope in-self calibration. This method does not need the support of external high-precision equipment compared with traditional calibration scheme and can be accomplished by user hand rotation. In this approach, Kalman filter is designed to perform the calibration procedure and estimate gyroscope bias error, scale factor error and non-orthogonal error. The system observability is analyzed and the dynamic rotating conditions under which the sensor errors become observable are derived. The design principles of optimal calibration procedure are provided as well. Both simulated and practical experiments are carried out to test the validation of the proposed calibration algorithm. The achieved results demonstrate that the introduced approach can provide promising calibration scheme for the low-cost MEMS gyroscope.


Introduction
Rapid development of MEMS (micro-electro-mechanical system) technology and advances in manufacturing industry has made it possible to produce low-cost consumer-grade inertial sensors (i.e., gyroscope, accelerometer). The MEMS inertial sensor is advantageous in chip-size minimization, low-cost manufacturing, lower power consumption, and has been applied in multiple applications, such as vehicle and pedestrian navigation, wearable electronic devices, augmented reality (AR). However, the MEMS inertial sensor suffers from various error sources (i.e., bias, scale factor error, non-orthogonal error, random noise) and thermal drifts, which can cause negative effects on its utilization. Hence, a calibration process is highly required to determine the sensor errors and mitigate the drift.
Calibration is defined as the process of comparing instrument outputs with known reference information to determine coefficients that force the output to agree with the reference information across the desired range of output values. A calibration experiment is conducted to compute deterministic errors of sensors in laboratory [1], and the calibrated parameters are employed to remove the sensor errors and derive a more reliable measurement. The calibration approaches and technologies have been well researched and studied, such as the local level frame (LLF) method and six-position static method [2]. However, these traditional calibration methods are primarily designed for in-lab tests and high-quality sensors, such as navigation or tactical grade IMUs (inertial measurement unit), and often require the use of special references (i.e., Earth rotation) such as alignment to a given frame or specialized equipment [3]. Therefore, using the specific high-level laboratory conditions to calibrate the low-cost MEMS sensor is always costly and meaningless.
Several researches or studies have been proposed to simplify the calibration procedure for MEMS inertial sensors. Skog [4] made use the fact that the norm of the gyroscope output ideally equals to the The main advantages and contributions of the proposed method include: (1) It is an easy, convenient and novel approach to calibrate 9 gyroscope error parameters during short term (approximately 3-5 min); (2) The system observability is analyzed, and the sensor error observable conditions are obtained. Hence, the optimal calibration experiment can be reasonably and theoretically designed.
(3) This approach does not require the support of external equipment and only needs to employ the pre-calibrated accelerometer output as reference signal. Moreover, it does not require other aiding information, such as GPS, magnetometer.

Main Error Sources of MEMS Gyroscope
Generally, gyroscope suffers from various error sources, such as bias, scale factor error, non-orthogonal error, acceleration errors (g-sensitivity), non-linearity, quantization error, and random errors [1]. Due to different manufacturing levels and chip production technologies, some kinds of sensor errors are dominant, while several others only have limited contribution. Therefore, it is critical to find out the main error sources of MEMS gyroscope and take them into consideration. We review the datasheet of different levels MEMS gyroscope (i.e., consumer level, high level), and summarize the sensor errors in Tables 2-4.  Table 2 lists the gyroscope characteristics of consumer-level (low level) IMU, MPU-9255/9150 [14,15]. Tables 3 and 4 provide the gyroscope characteristics of high-level inertial measurement units ADI16367 [16] and IMU-440 [17]. According to the comparison of different levels MEMS gyroscope performance, the bias, scale factor error and non-orthogonal error are considered as the main error sources. The bias is at least 1 • /s for the high-level MEMS gyroscope and is 5 • /s for the consumer-level sensor. An uncompensated gyroscope bias in the sensor measurement will introduce an angle error (in roll or pitch) that is proportional to time t. Accordingly, an error proportional to t 2 and t 3 will be involved in velocity and position respectively. The bias error compensation determines the solution performance and is an essential part in navigation algorithm.
The scale factor error is at least 1%, which means that it introduces 1 • /s error if the rotation velocity is 100 • /s and the error will enlarge with the increase of angular velocity. The non-orthogonal error, which denotes that the three axes are not perfectly aligned to the coordinate and the rotation in one axis will be deduced in the other two axes, is approximately 1 degree. In this assumption, if the three axes rotation velocities are 100 • /s, the non-orthogonal error will introduce 1.4 degrees measurement error to sensor output [3]. Hence, although the scale factor error and non-orthogonal error of high-level MEMS gyroscope is lower than that of consumer-level sensor, the measurement error caused by these errors is large and cannot be ignored.
For the nonlinearity error, it is 0.1% and is the least error among MEMS gyroscope error characteristic. The effect of the nonlinearity error is much less than that caused by other errors. For instance, its error is only 0.1 • /s when the rotation velocity is 100 • /s. Therefore, for the calibration algorithm proposed in this paper, we only consider the bias, scale factor error and non-orthogonal error, and can safely ignore the nonlinearity error.
In addition, this comparison also presents the effect of temperature on the gyroscope performance. For the consumer-level IMU, MPU-9255/9150, in the temperature range of −40 • C to +85 • C, the sensitivity scale factor variation is 4%, and the variation of Zero Rate Output (ZRO) is up to 20-30 • /s. The scale factor error change of ADI1367 over temperature is 40 ppm/ • C, and the bias variation with temperature is 0.01 • /s/ • C. The environmental conditions, such as the temperature variations and added voltage, will negatively impact the measurement performance of the MEMS gyroscope. Therefore, aiming to derive a more reliable navigation solution, it is more important to calibrate the MEMS gyroscope especially for the low-cost consumer-level, before usage in different application and environmental conditions.

Sensor Error Model
Gyroscopes are angular rate sensors that provide either angular rate or attitude depending on whether they are of the rate sensing or rate integrating type. Measurements of angular rate can be modeled by the observation equation as follows: where, ω b ib denotes gyroscope measurement vector. ω b ib is the true angular rate velocity vector. b g is the gyroscope instrument bias vector. S g is a matrix representing the gyroscope scale factor error. N g is a matrix representing non-orthogonal of the gyroscope triad. ε g is a vector representing the gyroscope sensor noise.
The matrices N g and S g are given as: where, S (·) are the scale factors for the three gyroscopes and G (·) represent the non-orthogonal error.

Bias Calibration
The bias is defined as the signal output of the sensor when it is not experiencing any rotation. According to this definition and the measurement error model in Equation (1), when the sensor is kept in stationary, the true angular rate velocity is zero and the items related to ω b ib in Equation (1) disappear. Thus, the gyroscope output is the bias error.
The random noise ε g involved in the gyroscope measurement, can be modeled as zero-mean white noise. The bias error can be calculated by averaging the gyroscope measurement in static period.
where, N denotes the collected measurement points. Generally, we can collect 3-5 min gyroscope measurement in static period and the bias error is calculated as the average.

Scale Factor and Non-Orthogonal Calibration
With given initial attitude derived from alignment procedure, the gyroscope measurement can be integrated to calculate the orientation information through Strapdown inertial navigation algorithm. However, the computed attitude will drift over time and the error is gradually accumulated because of the sensor error. In stationary or low dynamic condition, the accelerometer output can be used to estimate the orientation relative to horizontal plane (i.e., pitch and roll). The attitude derived from accelerometer output is independent in different time epochs and not affected by accumulated error. Hence, based on the different sensors' complimentary error propagation characteristics, we can make use of the accelerometer-derived attitude as reference signal to evaluate the attitude error introduced during integration process, and consequently determine the gyroscope error.
During the calibration process, the IMU is handheld by user and rotated along its axes slowly to avoid introducing external acceleration. The IMU orientation keeps varying during this procedure and the attitudes derived from different inertial sensors are compared to amend the attitude error and determine the sensor errors. A Kalman filter [18] is designed to estimate the scale factor and non-orthogonal errors of gyroscope. The attitude error propagation equation, which includes sensor error, is utilized as the system dynamic model. The relationship between the accelerometer output and attitude error is modeled as the measurement equation.

Dynamic Model
A simplified psi-angle error mode is applied as the dynamic model in Kalman filter, where the Psi-angle error represents the Euler angles between true navigation n frame (East-North-Up) and the platform frame (computed navigation frame n c frame). Adopting this kind of angle error mode is advantageous as it avoids non-linear problem and is easy to be implemented. The dynamic model is expressed as: where, ϕ denotes the attitude error expressed in n frame (i.e., navigation frame), ω n in denotes the n frame rotation angular rate vector relative to the inertial frame (i frame) expressed in n frame. C n b denotes the Direction Cosine Matrix (DCM) from body frame (i.e., b frame) to n frame. The symbol "×" denotes cross product of two vectors. δω b ib denotes the gyroscope output error. The items ω n in , δω n in in Equation (4) are: Because the hand-held rotation strategy is adopted to calibrate the sensor error, the velocity will be extremely low (approximately 0 m/s) and the position change is reasonable to be considered as zero. In addition, the velocity components and the position error are either divided by the Earth's radius (or its square) or are multiplied by the Earth rotation. Hence, the two rotation items are nearly zero, which is much less than the sensor error item, and as such can be safely ignored. Hence, the Equation (4) is simplified as: This equation illustrates that the attitude error is caused by angular velocity errors in measuring the rotation between the two frames (i.e., i frame and n frame). Due to the bias error has been calibrated by averaging the gyroscope output in static period, only the scale factor error and non-orthogonal error are considered and the equation is rewritten as: We model the scale factor error and non-orthogonal error in system dynamic model as random constant [19] and the dynamic model is described as: where, ϕ is 3 × 1 vector, which contains the pitch, roll and yaw angle errors δθ, δγ, δφ. O 3 × 3 denotes the 3 × 3 zero matrix. ω b ib × denotes the skew-symmetric matrix of gyroscope measurement. diag ω b ib denotes the diagonal matrix with the gyroscope measurement as its diagonal elements and w denotes the process noise and is set according to the gyroscope Angular Random Walk (ARW) [20].

Measurement Model
The acceleration residuals between the projection of accelerometer measurement in computed navigation frame n c and the local gravity acceleration are used to build the measurement model. The acceleration projection and gravity acceleration are written as: where, g n denotes the local gravity acceleration. a x , a y , a z denote the accelerometer measurements. denotes the misalignment between accelerometer and gyroscope traids. It should be noted that the accelerometer has been well calibrated and its biases, scale factor errors and non-orthogonal errors have been removed [21].
Followed by DCM chain rule, C n c b is expressed as: where, [ϕ×] denotes the skew matrix of attitude error. The acceleration residual is the difference between g n and f n c . Substitute Equation (10) into the residual and we can derive.
Therefore, the measurement model is derived as: where O 3 × 6 denotes a 3 × 6 zero matrix. H denotes the design matrix and v denotes the measurement noise.
The system model adopts the attitude error propagation equation, which combines the attitude error and sensor error. This model establishes a direct relationship between these two errors and is beneficial to increase the parameter observability. Moreover, the attitude error model is expressed in linear form and the Kalman filter can provide an optimal solution. The acceleration residuals, instead of the pitch and roll errors, are used as the observation vector. Although, using the orientation difference from various sensors is an intuitive way to implement the observation vector, it will introduce the singularity problem in the attitude calculation when the pitch angle reaches ± 90 • . However, the acceleration residuals include equivalent errors information of pitch and roll angles, without being affected by singularity problem.

Observability Definition
Observability describes the ability of estimating the system states [22], which is a measure for how well system states can be inferred by knowledge of the measurement. Analyzing the gyroscope sensor error observability is beneficial to know the state vector under which dynamic condition (i.e., continuous rotation, kept static) becomes observable and well estimated; thus, the gyroscope calibration scheme can be further optimally designed and implemented.  (14) where, A(t) and C(t) are, the n × n and p × n matrices whose entries are continuous functions of time defined over (−∞, +∞), respectively. The dynamic equation is observable at t 0 if there exists a finite time t 1 > t 0 such that for any state x 0 at time t 0 , the knowledge of the output y(t) over the time interval suffices to determine the state x 0 .

x(t) = A(t)x(t) y(t) = C(t)x(t)
Define a sequence of p × n observability matrix N 0 (t), N 1 (t), · · · N n−1 (t) by the equation: Suppose A(t) and C(t) in the system are analytic functions of t. Then the time-varying system A(t) and C(t) is observable at the time t 0 if there exists a finite time t 1 > t 0 such that the rank of the matrix is n.
The above observability tests are the same as finding a state vector x(t) such that For a time t 1 > t 0 , if there is no nonzero state that satisfies the above conditions, then the system is observable at t 0 .

Observability Analysis
According to observability definition, and the system transition and design matrices, F(t) and H, the system observability is determined by the rotational angular velocity ω and the IMU attitude (i.e., the transformation matrix C n b ). By substituting (9) and (13) into (15), the proposed system observability matrix is described as: where, O 2 × 3 denotes 2 × 3 zero matrix. C mn denotes the element in C n b matrix which is in the m th row and n th column. C n b denotes the transformation matrix between b frame and n frame and is written as: Sensors 2020, 20, 5430 9 of 20 The body frame (i.e.,b frame) in this paper is defined as Right-Forward-Up and the navigation frame (i.e., n frame) is defined as East-North-Up (ENU). Three consecutive rotations from the n frame to b frame is φ along z-axis, θ along x-axis and γ along x, y, z-axes. Rotating the IMU in hand in a free-style will lead that the system process equations become time-varying and complicate the theoretical observability analysis process. Aiming to derive a reliable analysis and estimation result, we consider the typical or special rotation and analyze the system observability in these conditions. The state vector observability is investigated under two dynamic conditions: (1) Rotation along vertical axis, which means that the IMU rotates along its sensitive axis when this axis is vertical to horizontal plan; (2) Rotation along horizontal axis, which means that the IMU rotates along its sensitive axis when it is aligning the horizontal plane.

Rotation along Vertical Axis
Assuming the IMU is put in flat platform with the x-axis and y-axis aligning to horizontal plane, and the z-axis is perpendicular to the plane. The rotation vector expressed in b frame is (0, 0, ω). The roll and pitch angles are both zero. In this condition, the items N 2 , N 3 . . . in observability matrix will become zero matrix and the observability matrix is written as: (17), and the observability test is described as: Correspondingly, only in the case that the elements δγ, δθ, G x , G y in state vector are zero, the condition listed in Equation (17) is satisfied. Hence, when the z-axis is vertical to the horizontal plane and the IMU is rotating with respect to it, the attitude error δγ, δθ and the non-orthogonal error G x , G y are observable.
Furthermore, when the y-axis is perpendicular to the horizontal plane and the IMU rotates along this axis, the attitude error and the non-orthogonal error G x , G z are observable. If the IMU is rotating with respect to x-axis and this axis is vertical to the horizontal plane, the attitude error and the non-orthogonal error G x , G z are observable.

Rotation along Horizontal Axis
Assuming the IMU is put in flat platform and the z-axis is perpendicular to the horizontal plane, the three axes rotation are separately (0, ω, 0). The initial roll and pitch angles are both zero. In this condition, we ignore the N 2 , N 3 . . . in observability matrix and the observability matrix is written as: Sensors 2020, 20, 5430 10 of 20 Consequently, the observability test is described as: The rank of observability matrix is 4, while five states are available in the observability test. The attitude error, δγ and δθ are not coupled with other elements, and these two states are always observable. The combination of the rest three items, y-axis scale factor S y and the non-orthogonal error G x , G z are observable. Because the non-orthogonal errors have become observable and can be well estimated when the IMU is rotating with respect to its vertical axis. Therefore, though the scale factor error is coupled with the non-orthogonal error, it still can be derived and estimated.
Accordingly, in the same condition, if the IMU rotates with respect to the x-axis, the combination of S y , G y , G z is observable. If the IMU has an initial 90 • or −90 • roll angle that its z-axis is aligned to the horizontal plane, the combination of S z , G x , G y is observable when the IMU rotates alongthe z-axis.

Observability Analysis Summarization
Based on the theoretical analysis of system observability, the relationship between the gyroscope error observability and dynamic conditions are summarized in the Table 5. Table 5 lists the observable states, orientations, and dynamic rotations. To summarize, the non-orthogonal error is observable when the IMU is rotating with respect to its vertical axis. The combination of three-sensor errors, including one scale factor error and two non-orthogonal errors, is observable when the IMU is rotating with respect to its horizontal axis.
The scale factor and non-orthogonal errors are coupled together when the IMU is rotating along the horizontal plane; however, the non-orthogonal error can be well estimated when the IMU rotates with respect to the vertical axis. Therefore, in order to acquire a reliable error estimation result, the calibration procedure should be designed as rotating along the vertical axis to derive the three non-orthogonal errors first, and then rotating with respect to the horizontal axis to estimate the three scale factor errors.

Observable States Pitch and Roll Angles Rotating Axis Rotation
Sensors 2020, 20, x FOR PEER REVIEW 11 of 21

Observable States Pitch and Roll Angles Rotating Axis Rotation
Combination

Test and Results
Simulated and practical experiments have been carried out to test and verify the proposed calibration algorithm. In simulated experiment, we first design the sensor rotation sequence to generate true sensor data according to observability analysis and add the preset sensor error to data which is used to perform the calibration approach. In practical experiment, the rotation follows the same sequence as that performed in simulation. In each experiment, the converge of estimated scale factor error and non-orthogonal error and their corresponding element in P matrix are drawn, compared and analyzed. In simulation experiment, we compare the estimated sensor error and preset error parameter to derive the absolute error and relative error. In practical experiment, the attitude calculated by sensor data with and without error compensation are compared to demonstrate the validation of calibrated parameter and all the figures are plotted in Matlab.

Test and Results
Simulated and practical experiments have been carried out to test and verify the proposed calibration algorithm. In simulated experiment, we first design the sensor rotation sequence to generate true sensor data according to observability analysis and add the preset sensor error to data which is used to perform the calibration approach. In practical experiment, the rotation follows the same sequence as that performed in simulation. In each experiment, the converge of estimated scale factor error and non-orthogonal error and their corresponding element in P matrix are drawn, compared and analyzed. In simulation experiment, we compare the estimated sensor error and preset error parameter to derive the absolute error and relative error. In practical experiment, the attitude calculated by sensor data with and without error compensation are compared to demonstrate the validation of calibrated parameter and all the figures are plotted in Matlab.

Test and Results
Simulated and practical experiments have been carried out to test and verify the proposed calibration algorithm. In simulated experiment, we first design the sensor rotation sequence to generate true sensor data according to observability analysis and add the preset sensor error to data which is used to perform the calibration approach. In practical experiment, the rotation follows the same sequence as that performed in simulation. In each experiment, the converge of estimated scale factor error and non-orthogonal error and their corresponding element in P matrix are drawn, compared and analyzed. In simulation experiment, we compare the estimated sensor error and preset error parameter to derive the absolute error and relative error. In practical experiment, the attitude calculated by sensor data with and without error compensation are compared to demonstrate the validation of calibrated parameter and all the figures are plotted in Matlab.

Test and Results
Simulated and practical experiments have been carried out to test and verify the proposed calibration algorithm. In simulated experiment, we first design the sensor rotation sequence to generate true sensor data according to observability analysis and add the preset sensor error to data which is used to perform the calibration approach. In practical experiment, the rotation follows the same sequence as that performed in simulation. In each experiment, the converge of estimated scale factor error and non-orthogonal error and their corresponding element in P matrix are drawn, compared and analyzed. In simulation experiment, we compare the estimated sensor error and preset error parameter to derive the absolute error and relative error. In practical experiment, the attitude calculated by sensor data with and without error compensation are compared to demonstrate the validation of calibrated parameter and all the figures are plotted in Matlab.

Simulation Experiment
Aiming to derive a reliable gyroscope error calibration result, the IMU is designed to rotate along vertical axis first to estimate the non-orthogonal error and consequently along horizontal axis to estimate the scale factor error based on the system observability analysis. The rotation scheme design is described in the Table 6.
The IMU keeps static with zero attitude in the first 10 s and rotates with respect to z-axis in clockwise and counterclockwise directions. Then the IMU rotates 90 degrees along y-axis to move the x-axis vertical to horizontal plane and rotates with respect to x-axis. Consequently, the IMU separately rotates along x, y, z-axes in horizontal plane. Table 6 lists the designed rotation sequence including the time period, direction and rotating axis.
The gyroscope error including scale factor error, non-orthogonal error, ARW, and the accelerometer error Velocity Random Walk (VRW) are simulated according to the consumer level MEMS sensor performance and listed in Table 7. These error parameters are substituted into the sensor measurement model to generate the output and are used for calibration.  Figure 1 shows the simulated inertial sensor measurement and Figure 1a,b separately illustrate the rotation and acceleration. The blue, red and yellow lines denote the measurement in x, y, z axes. The proposed calibration algorithm is employed to estimate the sensor error and the initial state vector and variance-covariance matrix P are set as:  The initial state is set as zero-vector due to the lack of sensor error prior knowledge. The P matrix is set based on the maximum sensor error. The error estimation result is drawn in Figure 2. The initial state is set as zero-vector due to the lack of sensor error prior knowledge. The P matrix is set based on the maximum sensor error. The error estimation result is drawn in Figure 2.  The initial state is set as zero-vector due to the lack of sensor error prior knowledge. The P matrix is set based on the maximum sensor error. The error estimation result is drawn in Figure 2.    Figure 2b shows the non-orthogonal error result. The blue, red and yellow line denote the error of x, y and z axes. Figure 3 shows the elements of P matrix, Figure 3a,b denote the elements corresponding to the scale factor error and the non-orthogonal error; Figure 3c,d denote the elements corresponding to attitude error and heading error. The reason for investigating the elements of covariance matrix P is that if its diagonal element experiences a large decrease from its initial value, the corresponding state becomes observable and its observability is large [20].
As shown in Figure 3, because the IMU keeps static in the first 10 s, the state vector is not observable, and the diagonal element of P matrix does not change. Starting from 10 s, the IMU rotates with respect to z-axis, the , converge to the reasonable value and their corresponding P matrix elements decrease dramatically. When the IMU rotates along x axis vertically beginning from 40 s, the begins to converge and its corresponding P matrix element decreases. The non-orthogonal   Figure 2b shows the non-orthogonal error result. The blue, red and yellow line denote the error of x, y and z axes. Figure 3 shows the elements of P matrix, Figure 3a,b denote the elements corresponding to the scale factor error and the non-orthogonal error; Figure 3c,d denote the elements corresponding to attitude error and heading error. The reason for investigating the elements of covariance matrix P is that if its diagonal element experiences a large decrease from its initial value, the corresponding state becomes observable and its observability is large [20].  For the scale factor error, the IMU rotates with respect to the x-axis and z-axis separately starting from 80 s and 126 s, then the and begin to converge and their corresponding P matrix elements experience a large decrease simultaneously. The scale factor error estimation result also demonstrates the observability analysis that rotating along the IMU horizontal axis is beneficial to make the scale factor error observable and estimate this error. As shown in Figure 3, because the IMU keeps static in the first 10 s, the state vector is not observable, and the diagonal element of P matrix does not change. Starting from 10 s, the IMU rotates with respect to z-axis, the G x , G y converge to the reasonable value and their corresponding P matrix elements decrease dramatically. When the IMU rotates along x axis vertically beginning from 40 s, the G z begins to converge and its corresponding P matrix element decreases. The non-orthogonal error estimation result verifies the observability analysis that rotating along vertical axis is able to make the non-orthogonal error observable and achieve a reliable calibration result.
For the scale factor error, the IMU rotates with respect to the x-axis and z-axis separately starting from 80 s and 126 s, then the S x and S z begin to converge and their corresponding P matrix elements experience a large decrease simultaneously. The scale factor error estimation result also demonstrates the observability analysis that rotating along the IMU horizontal axis is beneficial to make the scale factor error observable and estimate this error.
For the attitude error (i.e., pitch error, roll error), the two parameters are observable in the filter and is unrelated to the system dynamic condition. Hence, the P matrix elements corresponding to these two errors converge after the filter begins to work. The heading error is unobservable during the process, because the system measurement vector is uncorrelated to the heading error and cannot provide helpful correction, which explains that the its corresponding P matrix element diverges in the calibration procedure. In addition, the bias, scale factor error and non-orthogonal error estimation results are listed in Tables 8 and 9. The absolute error and relative error compared with true value are also provided. The gyroscope error calibration result is reasonable and reliable. The maximum relative error is 6.8%, and the minimum relative error is only 1%. Moreover, if the random noise of inertial sensor can be decreased in the simulated process, a more accurate result can be acquired.
Additionally, we utilize the singular value decomposition (SVD) of the observability stripped observation matrix to evaluate the degree of observability of each state. We separately calculate the singular value of the observability matrix with and without the observable state according to analysis result. Tables 10 and 11 give the result.  Tables 10 and 11 list the singular value of the SOM when the IMU rotates under different dynamic motion. In order to analyze the single state observability, we separately calculate the singular values with full states (9 states) and chosen states (6 or 7 states) and make the comparison. The zero singular values are ignored and only the non-zero values are listed in the tables.
As described in tables, two singular values 9.8 are available in each group. The reason is that the pitch error and roll error are always observable during the procedure, and 9.8 is their singular value. Hence, only need to analyze the other values. When the IMU is rotating with respect to the z-axis, the singular values of SOM (Stripped Observability Matrix) with the full states are 8.52 and 8.50. While the singular value become 0.61935 and 0.375294 when the states G x and G y are not considered in the observability matrix. It shows that the system observability degree decreases dramatically and only the attitude errors are observable in this condition. Meanwhile, it demonstrates that G x and G y are observable when the IMU rotates along the z-axis vertically. Based on this analysis and the comparison of singular values under various dynamic conditions, it comes to the same result with the observability analysis that the scale factor error is observable when the IMU rotates horizontally and the non-orthogonal error is observable when the IMU is rotating vertically.

Practical Experiment
We employ the consumer level MESM sensor MPU-9150 in practical experiment to test the proposed calibration algorithm. The MPU-9150 is a 9-axis inertial sensor which includes 3-axis gyroscope, 3-axis accelerometer and 3-axis magnetometer. The rotation sequence follows the designed scheme described in simulation experiment that rotating along the vertical axis first and then along the horizontal axis. The collected data is illustrated in Figure 4. gyroscope, 3-axis accelerometer and 3-axis magnetometer. The rotation sequence follows the designed scheme described in simulation experiment that rotating along the vertical axis first and then along the horizontal axis. The collected data is illustrated in Figure 4.
(a) Gyroscope measurement (b) Gyroscope measurement in static period   Figure 4 illustrates the collected inertial data in practical experiment. Figure 4a,c show the gyroscope and accelerometer measurement. Figure 4b shows the gyroscope measurement of static period in the first 10s. The rotation data of static period is averaged to derive the sensor bias. The initial state vector and covariance matrix P is set the same as that in the simulated experiment. The sensor error calibration result is shown in Figure 5.   Figure 4 illustrates the collected inertial data in practical experiment. Figure 4a,c show the gyroscope and accelerometer measurement. Figure 4b shows the gyroscope measurement of static period in the first 10s. The rotation data of static period is averaged to derive the sensor bias. The initial state vector and covariance matrix P is set the same as that in the simulated experiment. The sensor error calibration result is shown in Figure 5.   Figure 4 illustrates the collected inertial data in practical experiment. Figure 4a,c show the gyroscope and accelerometer measurement. Figure 4b shows the gyroscope measurement of static period in the first 10s. The rotation data of static period is averaged to derive the sensor bias. The initial state vector and covariance matrix P is set the same as that in the simulated experiment. The sensor error calibration result is shown in Figure 5.   Figure 5 shows the sensor error estimation result, Figure 5a denotes the scale factor error and Figure 5b shows the non-orthogonal error result. The blue, red and green lines denote the errors of x, y and z axes. Figure 6 shows the element of P matrix, Figure 6a,b denote the elements corresponding to the scale factor error and non-orthogonal error; Figure 6c,d denote the elements corresponding to the attitude error and heading error.
Due to the IMU keeps static in the first 10 s, the sensor errors (i.e., scale factor error, non-orthogonal error) are not observable and cannot be well estimated. After the static period, the IMU starts to rotate along z-axis vertically, the G x , G y begin to converge and their corresponding P matrix diagonal elements experience large decrease. Starting from 30 s, the IMU rotates along x-axis vertically and its corresponding P matrix decreases. The IMU starts to rotate along x, y, z axis separately approximately since 52 s, 72 s, 90 s, and the scale factor errors, S x , S y , S z also begin to gradually converge from these time epochs. Meanwhile, their corresponding P matrix diagonal elements decrease.
In order to test the validation of calibration algorithm, we rotate the IMU in free style and separately calculate the attitude using the sensor data with and without error compensation. The attitude is computed with the orientation integration algorithm and the results are compared with the reference. The attitude calculation results are illustrated in the following figures.
Sensors 2020, 20, x FOR PEER REVIEW 17 of 21 Figure 5 shows the sensor error estimation result, Figure 5a denotes the scale factor error and Figure 5b shows the non-orthogonal error result. The blue, red and green lines denote the errors of x, y and z axes. Figure 6 shows the element of P matrix, Figure 6a,b denote the elements corresponding to the scale factor error and non-orthogonal error; Figure 6c,d denote the elements corresponding to the attitude error and heading error.  Due to the IMU keeps static in the first 10 s, the sensor errors (i.e., scale factor error, non-orthogonal error) are not observable and cannot be well estimated. After the static period, the IMU starts to rotate along z-axis vertically, the , begin to converge and their corresponding P matrix diagonal elements experience large decrease. Starting from 30 s, the IMU rotates along x-axis vertically and its corresponding P matrix decreases. The IMU starts to rotate along x, y, z axis separately approximately since 52 s, 72 s, 90 s, and the scale factor errors, , , also begin to gradually converge from these time epochs. Meanwhile, their corresponding P matrix diagonal elements decrease.
In order to test the validation of calibration algorithm, we rotate the IMU in free style and separately calculate the attitude using the sensor data with and without error compensation. The attitude is computed with the orientation integration algorithm and the results are compared with the reference. The attitude calculation results are illustrated in the following figures. Figure 7 shows the attitude calculation result, Figure 7a,c,e show the roll, pitch, heading results and Figure 7b,d,f show the attitude results in the selected periods for better comparison. The blue line denotes the attitude result calculated by the inertial sensor data which has been compensated by the error parameters derived from proposed calibration approach. The red line denotes the attitude calculated by the data in which only the bias has been removed, and the green line denotes the reference.   Figure 7b,d,f show the attitude results in the selected periods for better comparison. The blue line denotes the attitude result calculated by the inertial sensor data which has been compensated by the error parameters derived from proposed calibration approach. The red line denotes the attitude calculated by the data in which only the bias has been removed, and the green line denotes the reference.   As shown in the right side three figures which have been enlarged in selected periods, the attitude with error compensation (blue line) is closer to the reference (green line) and the attitude without error compensation (red line) drifts away from the reference. The attitude comparison demonstrates the validation of gyroscope error estimation result because the attitude performance has been improved through well sensor error compensation, while the attitude without error compensation has a poorer performance. The attitude error result and its statistical result are shown in Figure 8 and Table 12  As shown in the right side three figures which have been enlarged in selected periods, the attitude with error compensation (blue line) is closer to the reference (green line) and the attitude without error compensation (red line) drifts away from the reference. The attitude comparison demonstrates the validation of gyroscope error estimation result because the attitude performance has been improved through well sensor error compensation, while the attitude without error compensation has a poorer performance. The attitude error result and its statistical result are shown in Figure 8 and Table 12.  Figure 8 shows the attitude error result, Figure 8a-c separately denote the pitch, roll and heading errors. The blue and red line respectively denote the attitude calculated with and without sensor error compensation. Table 12 lists the error statistical result of the two groups attitude result and illustrates that the attitude result with error compensation has an overall better performance.
As shown in Figure 8, in the first 8 s during which the IMU keeps static, the attitude errors of the two groups are approximately same and almost equal to zero. Then when the IMU begins to rotate, the attitude error curves become different and it is obvious that the attitude with error compensation performed well than that without error compensation.
During the static period, due to the bias error has been well estimated and removed according to the gyroscope output average, the performances of two groups attitude are nearly same. As analyzed in the MEMS sensor error source section, the scale factor and non-orthogonal errors are proportional to the angular velocity. Therefore, during the rotating stage, these two errors begin to affect the calculation result, and the inertial data with error compensation constrains the negative effect and leads to a more reliable solution. While the sensor data without error compensation suffers from the error effect and the attitude calculated by this group data will drift faster and has a poor performance.   Figure 8 shows the attitude error result, Figure 8a-c separately denote the pitch, roll and heading errors. The blue and red line respectively denote the attitude calculated with and without sensor error compensation. Table 12 lists the error statistical result of the two groups attitude result and illustrates that the attitude result with error compensation has an overall better performance.
As shown in Figure 8, in the first 8 s during which the IMU keeps static, the attitude errors of the two groups are approximately same and almost equal to zero. Then when the IMU begins to rotate, the attitude error curves become different and it is obvious that the attitude with error compensation performed well than that without error compensation.
During the static period, due to the bias error has been well estimated and removed according to the gyroscope output average, the performances of two groups attitude are nearly same. As analyzed in the MEMS sensor error source section, the scale factor and non-orthogonal errors are proportional to the angular velocity. Therefore, during the rotating stage, these two errors begin to affect the calculation result, and the inertial data with error compensation constrains the negative effect and leads to a more reliable solution. While the sensor data without error compensation suffers from the error effect and the attitude calculated by this group data will drift faster and has a poor performance.
The proposed calibration algorithm is implemented through Kalman filter. The dynamic and measurement model are all expressed in linear form, the state and measurement vector dimension The proposed calibration algorithm is implemented through Kalman filter. The dynamic and measurement model are all expressed in linear form, the state and measurement vector dimension are separately 9 and 3 which is not large; thus, the approach's complexity is O(9) and does not occupy much computation resource. In addition, sequential measurement update can be applied in filter to avoid calculating the inverse of matrix. Hence, the algorithm only needs limited computing resources and can be embedded in sensor's internal processor to implement error calibration without support of external CPU.

Conclusions
This paper presents a real-time low-cost MEMS gyroscope calibration method without external equipment. This calibration algorithm is able to identify the sensor bias, scale factor error and non-orthogonal error, and the calibration experiment can be easily accomplished by user hand-hold rotation. The system observability is analyzed, and the observable condition of sensor error is investigated. This research can guide the calibration scheme design to save time and achieve reliable estimation result. Both simulated and practical experiments have demonstrated the validation of proposed calibration approach.
Author Contributions: Q.Z. conceived the idea, designed the Kalman filter in algorithm, implemented the proposed software for verification, performed the experiment and wrote this paper. G.Y. supervised this research work and provided much valuable revising advice. H.L. and N.Z. helped to carry out experiments, review this paper and proposed some revising suggestions. All authors have read and agreed to the published version of the manuscript.