Complete Tri-Axis Magnetometer Calibration with a Gyro Auxiliary

Magnetometers combined with inertial sensors are widely used for orientation estimation, and calibrations are necessary to achieve high accuracy. This paper presents a complete tri-axis magnetometer calibration algorithm with a gyro auxiliary. The magnetic distortions and sensor errors, including the misalignment error between the magnetometer and assembled platform, are compensated after calibration. With the gyro auxiliary, the magnetometer linear interpolation outputs are calculated, and the error parameters are evaluated under linear operations of magnetometer interpolation outputs. The simulation and experiment are performed to illustrate the efficiency of the algorithm. After calibration, the heading errors calculated by magnetometers are reduced to 0.5° (1σ). This calibration algorithm can also be applied to tri-axis accelerometers whose error model is similar to tri-axis magnetometers.


Introduction
Strap-down inertial measurement units (SIMU) which consist of gyroscopes and accelerometers, have been applied to orientation estimation for a long time [1]. Nowadays, magnetometers rigidly mounted to SIMU, also known as MARG (Magnetic, Angular Rate, and Gravity), are popular for achieving better performance. However, the magnetometer outputs suffer from magnetic field distortions and sensor errors, including the hard iron effect, soft iron effect, zero bias error, scale factor error, non-orthogonal error, misalignment error, and noise [2]. Magnetometer calibrations are always advisable to be done prior to application.
The calibration methods which generally rely on non-magnetic turntables [3][4][5] and Helmholtz coils [6][7][8][9], are proposed first. The magnetometer error parameters can be precisely calibrated. However, the cumbersome processes and calibration platforms bring about the issues of inconvenience and expense. Many stand-alone magnetometer calibration algorithms based on an ellipse-fitting method have been put forward [10][11][12][13][14][15]. These algorithms are low-cost and easy to use. The aim of these ellipse-fitting based algorithms is to compensate the magnetometer output data from lying on an ellipsoid to a sphere. However, the misalignment error between magnetometers and assembled platforms, which causes a rotation of the sphere, is unable to be calibrated. For complete tri-axis magnetometer calibration, a series of calibration algorithms in combination with accelerometers [16][17][18][19][20][21] are proposed. First, the ellipsoid-fitting algorithms are performed to calibrate errors, except the misalignment error. Then, optimization problems based on the inner products of the tri-axis magnetometer outputs and accelerometer outputs are established and solved for misalignment error estimation. The author further applies gyroscopes with magnetometers and accelerometers to run an EKF (Extended Kalman Filter) for a maximum likelihood calibration process [22,23]. The nonlinear and iterative operations as well as the reliance on initial values create the drawback of heavy computation burdens for these algorithms. Tri-axis magnetometers and gyroscopes are also simultaneously applied in [24,25] to establish an EKF and an adaptive identification algorithm for magnetometer calibration without the use of an accelerometer. The error models of a tri-axis magnetometer are simplified. Only parts of the errors are calibrated, which degrade the calibration performance.
In this work, a tri-axis gyroscope is applied as an auxiliary for complete tri-axis magnetometer calibration. Firstly, the sensor set of the tri-axis magnetometer and the tri-axis gyroscope is mounted into a non-magnetic cuboid frame. Then, the frame is put on a table and rotated for at least one cycle for each side of the frame, respectively. Finally, the complete tri-axis magnetometer calibration is achieved through linear algebra operations of the magnetometer and gyroscope output data. After calibration, the heading errors calculated by magnetometers are reduced to 0.5 • (1σ). This algorithm requires a cuboid frame to mount the sensor set, which is easy to obtain and inexpensive. The proposed method, which employs the linear algebra operations, is of low computational complexity. Furthermore, because the error model of tri-axis accelerometers has a similar form to tri-axis magnetometers, this calibration algorithm can also be applied to tri-axis accelerometer calibrations.

Error Model
To illustrate the calibration algorithm, the error model of a tri-axis magnetometer is established. Similar models have been established in [10][11][12][13][14][15][16]. For convenience, the complete tri-axis error model is discussed again with misalignment errors included. The hysteresis errors, temperature-dependent errors, and time-varying errors will not be discussed in this paper. The output errors are classified into two categories: magnetic field distortions and sensor errors.

Magnetic Field Distortions
The magnetic field distortions are due to the ferromagnetic materials attached to the sensor frame, including the hard iron effect and soft iron effect. The hard iron effect is a constant additional magnetic field produced by permanent magnets on the assembled platform. It is represented as a 3 × 1 vector, denoted as b H . The soft iron effect is produced by the magnetization of soft magnets. It exerts effects on the magnitude and orientation of the magnetic field according to the external field orientation, which is represented as a 3 × 3 matrix, denoted as C si . The i-row, j-column element in the matrix indicates the influence of the external j-direction field to the i-direction field.

Sensor Errors
The sensor errors of the tri-axis magnetometer include the zero bias error, scale factor error, non-orthogonal error, misalignment error, and noise. These errors are mainly due to machining and installation defects. The zero bias error makes constant offsets to each axis of the magnetometer, which is represented as a 3 × 1 vector, denoted as b zb . The scale factor error comes from the sensitivity inconsistencies of each sensor. It is represented as a 3 × 3 diagonal matrix, denoted as C sf . The elements on the diagonal represent the sensitivities of each axis sensor. The non-orthogonality between individual sensors introduces an inter-axis coupling output error, which is represented as a 3 × 3 upper triangular matrix, denoted as C no . In general conditions of small angle errors, the diagonal elements of C no are close to 1. The misalignment error represents the angular misalignment between the tri-axis magnetometer set and the assembled platform, which is represented as a 3 × 3 unit rotation matrix, denoted as C m . It can be simplified into an anti-symmetric matrix, with diagonal elements equal to 1. The output noise of individual sensor is assumed as uncorrelated Gaussian white noise. The output noise is represented as a 3 × 1 vector, denoted as ε.
Above all, the error model of the tri-axis magnetometer is given as: In the equation, m stands for the magnetometer output vector. K stands for the total error coefficient matrix, and b stands for the total bias error vector. Note that in this paper,~stands for ideal quantities,ˆstands for interpolation quantities, and stands for compensated quantities.

Calibration Algorithm Development
The calibration algorithm is developed for evaluating the error coefficient matrix K and the bias error vector b. In this calibration algorithm, the sensor set is firstly mounted into a non-magnetic cuboid frame. The six sides of the cuboid are numbered from 1 to 6. Side 1 is opposite to side 2. Side 3 is opposite to side 4, and side 5 is opposite to side 6. The sensor set is aligned with the frame. The z axis of the sensor set is perpendicular to side 1 and side 2. The y axis of the sensor set is perpendicular to side 3 and side 4. The x axis of the sensor set is perpendicular to side 5 and side 6, as shown in Figure 1. The calibration process is divided into two steps. In the first step, the tri-axis magnetometer interpolation outputs at specified spin angular intervals are calculated. In the second step, the error coefficient matrix and bias error vector are calculated through linear algebra operations of magnetometer interpolation outputs.
In the equation, m stands for the magnetometer output vector. K stands for the total error coefficient matrix, and b stands for the total bias error vector. Note that in this paper, ~ stands for ideal quantities, ^ stands for interpolation quantities, and ͡ stands for compensated quantities.

Calibration Algorithm Development
The calibration algorithm is developed for evaluating the error coefficient matrix K and the bias error vector b. In this calibration algorithm, the sensor set is firstly mounted into a non-magnetic cuboid frame. The six sides of the cuboid are numbered from 1 to 6. Side 1 is opposite to side 2. Side 3 is opposite to side 4, and side 5 is opposite to side 6. The sensor set is aligned with the frame. The z axis of the sensor set is perpendicular to side 1 and side 2. The y axis of the sensor set is perpendicular to side 3 and side 4. The x axis of the sensor set is perpendicular to side 5 and side 6, as shown in Figure 1. The calibration process is divided into two steps. In the first step, the tri-axis magnetometer interpolation outputs at specified spin angular intervals are calculated. In the second step, the error coefficient matrix and bias error vector are calculated through linear algebra operations of magnetometer interpolation outputs.

Magnetometer Interpolation Output Calculation
The cuboid frame is put on a table and rotated over one cycle by hand. This is repeated for each numbered side facing upwards. Then the tri-axis magnetometer output data m k and gyroscope output data ω k are obtained, k = 1-6, where k stands for the number of the upward side. The process is shown in Figure 3. Here are some points to note: (1) the sequence of which side is facing up can be arbitrary; (2) the orientation of the flank when putting the frame on the table can be arbitrary (finding north is not required); (3) both clockwise and anticlockwise rotations are accepted; (4) constant rotation speed is not required; and (5) the number of turns must be larger than 1. The tri-axis magnetometer interpolation outputs at specified spin angular intervals are calculated from m k and ω k . The interpolation flowchart is shown in Figure 3. Firstly, the number of turns n, and the angular interval Δθ are set. The value of 360/Δθ must be an integer. The upward side number k is to be verified. Next, by integrating the output data ω ⊥ k of the gyroscope whose sensitive

Magnetometer Interpolation Output Calculation
The cuboid frame is put on a table and rotated over one cycle by hand. This is repeated for each numbered side facing upwards. Then the tri-axis magnetometer output data m k and gyroscope output data ω k are obtained, k = 1-6, where k stands for the number of the upward side. The process is shown in Figure 2. Here are some points to note: (1) the sequence of which side is facing up can be arbitrary; (2) the orientation of the flank when putting the frame on the table can be arbitrary (finding north is not required); (3) both clockwise and anticlockwise rotations are accepted; (4) constant rotation speed is not required; and (5) the number of turns must be larger than 1.
In the equation, m stands for the magnetometer output vector. K stands for the total error coefficient matrix, and b stands for the total bias error vector. Note that in this paper, ~ stands for ideal quantities, ^ stands for interpolation quantities, and ͡ stands for compensated quantities.

Calibration Algorithm Development
The calibration algorithm is developed for evaluating the error coefficient matrix K and the bias error vector b. In this calibration algorithm, the sensor set is firstly mounted into a non-magnetic cuboid frame. The six sides of the cuboid are numbered from 1 to 6. Side 1 is opposite to side 2. Side 3 is opposite to side 4, and side 5 is opposite to side 6. The sensor set is aligned with the frame. The z axis of the sensor set is perpendicular to side 1 and side 2. The y axis of the sensor set is perpendicular to side 3 and side 4. The x axis of the sensor set is perpendicular to side 5 and side 6, as shown in Figure 1. The calibration process is divided into two steps. In the first step, the tri-axis magnetometer interpolation outputs at specified spin angular intervals are calculated. In the second step, the error coefficient matrix and bias error vector are calculated through linear algebra operations of magnetometer interpolation outputs.

Magnetometer Interpolation Output Calculation
The cuboid frame is put on a table and rotated over one cycle by hand. This is repeated for each numbered side facing upwards. Then the tri-axis magnetometer output data m k and gyroscope output data ω k are obtained, k = 1-6, where k stands for the number of the upward side. The process is shown in Figure 3. Here are some points to note: (1) the sequence of which side is facing up can be arbitrary; (2) the orientation of the flank when putting the frame on the table can be arbitrary (finding north is not required); (3) both clockwise and anticlockwise rotations are accepted; (4) constant rotation speed is not required; and (5) the number of turns must be larger than 1. The tri-axis magnetometer interpolation outputs at specified spin angular intervals are calculated from m k and ω k . The interpolation flowchart is shown in Figure 3. Firstly, the number of turns n, and the angular interval Δθ are set. The value of 360/Δθ must be an integer. The upward side number k is to be verified. Next, by integrating the output data ω ⊥ k of the gyroscope whose sensitive The tri-axis magnetometer interpolation outputs at specified spin angular intervals are calculated from m k and ω k . The interpolation flowchart is shown in Figure 3. Firstly, the number of turns n, and the angular interval ∆θ are set. The value of 360/∆θ must be an integer. The upward side number k is to be verified. Next, by integrating the output data ω k ⊥ of the gyroscope whose sensitive axis is perpendicular to the table, the rotation angles θ k of each sample points are derived. After which, the interpolation outputsm k are derived from m k and θ k . Then, the upward side k is changed to a different number, and the program is repeated. Finally, the interpolation outputsm k i are obtained, where i = 1, 2, . . . , 360n/∆θ, k = 1, 2, . . . , 6. axis is perpendicular to the table, the rotation angles θ k of each sample points are derived. After which, the interpolation outputs m k are derived from m k and θ k . Then, the upward side k is changed to a different number, and the program is repeated. Finally, the interpolation outputs m i k are obtained, where i = 1, 2,…, 360n/Δθ, k = 1, 2,…, 6.

Error Coefficient Calculation
The geomagnetic field is stable. It can hardly change within days and is assumed to be constant in the calibration process. The geomagnetic field component perpendicular to the table is denoted as h ⊥ . When side 1 faces upwards, m z 1 is equal to h ⊥ , and Equation (1) can be rewritten as: where C i z is the ideal rotation matrix around z axis of sensor set, which is expanded as: cos sin 0 Considering that: by substituting Equation (4) into Equation (2), there is:

Error Coefficient Calculation
The geomagnetic field is stable. It can hardly change within days and is assumed to be constant in the calibration process. The geomagnetic field component perpendicular to the table is denoted as h ⊥ . When side 1 faces upwards, ∼ m 1 z is equal to h ⊥ , and Equation (1) can be rewritten as: where C z i is the ideal rotation matrix around z axis of sensor set, which is expanded as: Considering that: by substituting Equation (4) into Equation (2), there is: Similarly, when side 2 faces upwards, When side 3 faces upwards, ∼ m 3 y is equals to h ⊥ . There are: When side 4 faces upwards, ∼ m 4 y is equals to −h ⊥ . There are: When side 5 faces upwards, ∼ m 5 x is equals to h ⊥ . There are: When side 6 faces upwards, ∼ m 6 x is equals to −h ⊥ . There are: By Equations (2)-(15), the error coefficient matrix K and bias error vector b are solved through linear operations as follows: 1.
By linear operation of Equations (8)- (11), h ⊥ K 2 is solved: 4. By linear operation of Equations (12)-(15), h ⊥ K 1 is solved: From Equations (16)- (19), h ⊥ K and b are derived. h ⊥ is a constant, which makes a scale for the tri-axis outputs. This scale has no effect on the orientation estimations. In many applications, the magnetometer outputs are normalized. Therefore, the tri-axis outputs can be compensated as: However, in some applications, such as magnetic field measurements, h ⊥ needs to be determined. h ⊥ can be obtained by geomagnetic field models, or by a single-axis calibrated magnetometer. In general cases, the non-diagonal elements of K are close to zero, so h ⊥ can also be evaluated as: Finally, the tri-axis magnetometer outputs are compensated as follows: The program code of the calibration algorithm is presented in the Appendix A.

Calibration Error Analysis
The gyro errors, cuboid frame errors, and rotation operation errors may degrade the calibration performance. The calibration errors introduced by the gyro bias, gyro misalignment, non-orthogonality of the cuboid frame, and rotation deviation are analyzed, respectively.

Gyro Bias
Gyro bias is a main error of gyroscopes. Most of the bias error can be estimated during stationary time periods, and the residual tri-axis gyroscope bias error is assumed to be: When side 1 faces upwards, the magnetometer interpolation outputs are: The actual rotation matrix around the z axis is: ∆t is the time interval between two adjacent interpolation points. The angle measurements by integrating the gyroscope output data are of good accuracy in the short-term, and the angle errors are of small value. The total time during one direction rotation is assumed to be less than 10 s, and by testing the gyroscope in the experiment, the bias error is less than 0.5 • /s, which means the angle error estimated by the gyroscope is less than 5 • . Under the condition of small angle errors, Equation (25) is rearranged as follows: where ∆C z i is the rotation matrix error caused by the gyro bias at i-th interpolation point, given as: Λ 1 is defined as the sum of the rotation matrix errors caused by the gyro bias when side 1 faces upwards, derived as: Similarly, there is: Thus, the calibration error of b caused by gyro bias is derived by applying Equation (29) into Equation (16), as shown in Equation (30): If the initial orientations of the flank when side k and side k + 1 (k = 1, 3, 5) face upwards are the same, the gyro bias error would cause little calibration error to b, and if the initial orientation of the flank are opposite, the gyro bias error would cause maximal calibration error to b.
The calibration error of K caused by the gyro bias is derived by applying Equation (29) into Equations (17)- (19), as shown in Equation (31): If the initial orientations of the flank when side k and side k + 1 (k = 1, 3, 5) face upwards are opposite, the gyro bias error would cause little calibration error to K, and if the initial orientation of the flank are the same, the gyro bias error would cause maximal calibration error to K.
∆b gb and ∆K gb satisfy: where F is the total geomagnetic field intensity, and I is the geomagnetic inclination angle. Considering I = 59 • , ∆t = 0.02 s, ∆θ = 1 • , b g = 0.5 • /s, and that the magnetometer outputs are normalized, even under the worst initial orientation conditions, the elements of ∆b gb and ∆K gb are less than 0.0034 and 0.0060, respectively.

Gyro Misalignment
When installing the gyroscope into the cuboid frame, the misalignment between gyro axes and cuboid faces could be introduced. This would affect the gyro outputs, and cause estimation errors of the rotation angle. The gyro misalignment is presented as a rotation matrix of the small angle: o(·) denotes the infinitesimal. The gyro output errors are: ∆g 1 , ∆g 2 are the gyro output errors whose axes are parallel to the table, and ∆g ⊥ is the gyro output error whose axis is perpendicular to the table. The misalignment angle causes little error on g ⊥ , which is applied to rotation estimation.
Thus, the gyro misalignment error introduces slight calibration errors, which can be ignored.

Non-Orthogonality of Cuboid Frame
Non-orthogonality of the cuboid frame is indicated in Figure 4. The α, β and γ are the deviation angles, which are of small value.
o(•) denotes the infinitesimal. The gyro output errors are: Δg 1 , Δg 2 are the gyro output errors whose axes are parallel to the table, and Δg ⊥ is the gyro output error whose axis is perpendicular to the table. The misalignment angle causes little error on g ⊥ , which is applied to rotation estimation.
Thus, the gyro misalignment error introduces slight calibration errors, which can be ignored.

Non-Orthogonality of Cuboid Frame
Non-orthogonality of the cuboid frame is indicated in Figure 4. The α, β and γ are the deviation angles, which are of small value. When sides 1-6 are perpendicular to the table, respectively, the output errors of the gyroscope whose axis is perpendicular to the table are derived, respectively: It can be deduced from Equation (35) that the non-orthogonality of the cuboid frame has little effect on the vertical gyro outputs.
The magnetometer output error, caused by the non-orthogonality of the cuboid frame, is represented when side 1 faces upwards: In the same way, the magnetometer output errors with each side upwards are derived as: When sides 1-6 are perpendicular to the table, respectively, the output errors of the gyroscope whose axis is perpendicular to the table are derived, respectively: It can be deduced from Equation (35) that the non-orthogonality of the cuboid frame has little effect on the vertical gyro outputs.
The magnetometer output error, caused by the non-orthogonality of the cuboid frame, is represented when side 1 faces upwards: In the same way, the magnetometer output errors with each side upwards are derived as: Thus, by applying Equations (37) into Equation (16), the calibration error of b caused by the non-orthogonality of the cuboid frame is calculated: By applying Equation (37) into Equations (17)- (19), the calibration error of K is calculated: The non-orthogonality error of the cuboid frame introduces calibration errors to four elements of K, and does not affect the other elements of K and b. The errors are equal to the deviation angles. Assuming the machining verticality is 0.02, which is easy to be satisfied, and the side length of the cuboid frame is 100 mm, the errors to the four elements of K are less than 2 × 10 −4 , which could be ignored.

Rotation Deviation
While rotating the frame, the orientation of the rotations may change slightly due to the table being uneven and vibrations caused by friction during the rotation. This can result in unwanted effects on the gyro outputs and magnetometer outputs. The orientation deviation is presented as: The gyro output error caused by the rotation deviation is: g ⊥ is slightly affected by the rotation deviation. When side 1 faces upwards, the magnetometer output errors caused by the rotation deviation is: In the same way, the magnetometer output errors when the other sides face upwards are obtained: By applying Equation (43) into Equation (16), the calibration error of b caused by the rotation deviation is calculated: By applying Equation (43) into Equations (17)- (19), the calibration error of K is calculated: The deviation angles ε k j,i are assumed to be normally distributed with zero mean and variance of Q ε , and the ε k j,i and m k j,i are uncorrelated. The variance of the calibration error caused by the rotation deviation satisfies: The magnetometer outputs are considered as being normalized, such that the norm of the magnetic field F is set as 1. If √ Q = 1 • , I = 59 • , ∆θ = 1 • , and n = 2, the standard deviations of the elements in ∆b rd are less than 1.2 × 10 −3 , and the standard deviations of diagonal elements and off-diagonal elements in ∆K rd are less than 2.2 × 10 −3 and 2.6 × 10 −3 , respectively.
Above all, the calibration errors introduced by the gyro bias, gyro misalignment, non-orthogonality of the cuboid frame, and rotation deviation are slight and acceptable under normal conditions.

Simulation
To illustrate the calibration accuracy and computational requirements of the proposed method, the ellipsoid-fitting calibration method in [13] is performed for comparison. Hereafter, the ellipsoid-fitting calibration method in [13] is referred to as the EM, to distinguish from the proposed calibration method is referred to as the PM.
The local magnetic field is acquired from the world magnetic model [26]. As an example, the north, east, and down components of the Beijing area are 27,959.6 nT, −3296.3 nT, and 46,697.7 nT, respectively. Considering the sensor set used in the following experiment, the magnetometer output noise is set to 200 nT (1σ). The gyroscope output noise is 0.05 • /s (1σ) and the zero drift is 0.5 • /s for each individual sensor. The sample frequency is 100 Hz. The number of turns n is set to 2, and the interpolation angular interval ∆θ is set to 1 • .
Two error conditions of the tri-axis magnetometer are considered for further comparison. In the error condition 1, all of the tri-axis magnetometer errors are included. In the error condition 2, the misalignment error, soft iron error, and non-orthogonality error of the magnetometer, which cause rotations of the output data and degrade the calibration performance of the EM, are removed.

Error Condition 1
The magnetometer errors are set and shown in Table 1. Table 1. Tri-axis magnetometer errors.

Error
Value Error Value The tri-axis magnetometer output data during the calibration process are plotted in Figure 5, where RD refers to raw data, and ID refers to ideal data. After calibration using the PM, the calibrated outputs are consistent with the ideal outputs. The distorted ellipses are rectified into circles, centered at [0, 0], and parallel to the sides of the cuboid, which means the hard iron effect, soft iron effect, zero bias error, scale factor error, non-orthogonal error, and misalignment error are compensated. After calibration using the EM, as shown in Figure 6, the distorted ellipses are rectified into circles, and the calibrated outputs are distributed on the surface of a sphere centered at [0, 0, 0], which means that the total scale error and total zero bias error, introduced by the scale factor error, zero bias error, soft iron error, and non-orthogonality error, are compensated. However, the calibrated outputs have a rotation error from the ideal outputs. The rotation error has the same matrix form as the misalignment error, which is considered as a total misalignment error introduced by the misalignment error, soft iron error, and non-orthogonality error [12]. Using the arctangent function, the heading angles are calculated from the output data of twoaxis magnetometers parallel to the table. When side 1 faces upwards, the heading angle is calculated by the x-and y-axis magnetometer outputs. When side 3 faces upwards, the heading angle is calculated by the x-and z-axis magnetometer outputs. When side 5 faces upwards, the heading angle is calculated by the y-and z-axis magnetometer outputs. By comparing the raw and calibrated angle errors under the conditions of side 1, side 3, and side 5 facing upwards, respectively, the calibration performances of the EM and PM are illustrated, as shown in Figure 7. Using the arctangent function, the heading angles are calculated from the output data of two-axis magnetometers parallel to the table. When side 1 faces upwards, the heading angle is calculated by the xand y-axis magnetometer outputs. When side 3 faces upwards, the heading angle is calculated by the xand z-axis magnetometer outputs. When side 5 faces upwards, the heading angle is calculated by the yand z-axis magnetometer outputs. By comparing the raw and calibrated angle errors under the conditions of side 1, side 3, and side 5 facing upwards, respectively, the calibration performances of the EM and PM are illustrated, as shown in Figure 7.  After calibration using the EM, the angle errors are not reduced, but such errors become more consistent when different sides face upwards. After calibration using the PM, the angle errors are limited to 0.5° (1σ), with the maximum less than 1.5°, as shown in Table 2.  After calibration using the EM, the angle errors are not reduced, but such errors become more consistent when different sides face upwards. After calibration using the PM, the angle errors are limited to 0.5° (1σ), with the maximum less than 1.5°, as shown in Table 2. After calibration using the EM, the angle errors are not reduced, but such errors become more consistent when different sides face upwards. After calibration using the PM, the angle errors are limited to 0.5 • (1σ), with the maximum less than 1.5 • , as shown in Table 2.

Error Condition 2
Under error condition 2, the misalignment error C m , soft iron error C si , and non-orthogonality error C no of the magnetometer are removed and the other parameters keep the same. The simulation program is repeated under this condition. The angle errors are eliminated effectively both using the EM and the PM, as shown in Figure 8 and Table 3. When there is no misalignment error, soft iron error, and non-orthogonality error, this shows that the EM performs efficiently. The PM performs efficiently under both error conditions.

Error Condition 2
Under error condition 2, the misalignment error Cm, soft iron error Csi, and non-orthogonality error Cno of the magnetometer are removed and the other parameters keep the same. The simulation program is repeated under this condition. The angle errors are eliminated effectively both using the EM and the PM, as shown in Figure 8 and Table 3. When there is no misalignment error, soft iron error, and non-orthogonality error, this shows that the EM performs efficiently. The PM performs efficiently under both error conditions.  The amounts of computation of the EM and PM are compared. Both the EM and PM are divided into two steps. In the EM, step 1 is to perform ellipsoid-fitting using least squares technique, and step 2 is to extract the error coefficients from the ellipsoid parameters. In the PM, step 1 is to calculate the magnetometer interpolation outputs using the linear interpolation method, and step 2 is to calculate the error coefficients from the interpolation outputs by linear operations. The number of data for ellipsoid fitting is set to be equal to the number of interpolation outputs, which is equal to 360n/Δθ.  The amounts of computation of the EM and PM are compared. Both the EM and PM are divided into two steps. In the EM, step 1 is to perform ellipsoid-fitting using least squares technique, and step 2 is to extract the error coefficients from the ellipsoid parameters. In the PM, step 1 is to calculate the magnetometer interpolation outputs using the linear interpolation method, and step 2 is to calculate the error coefficients from the interpolation outputs by linear operations. The number of data for ellipsoid fitting is set to be equal to the number of interpolation outputs, which is equal to 360n/∆θ. The simulation program is repeated ten times under both of the error conditions using Matlab, and the average computation time of these two algorithms is shown in Table 4. The step 1 in the EM requires large amounts of computation. The linear operations in the PM requires small computational amounts. The EM requires computations amounting to more than 10 times that of the PM. Table 4. Computation time (simulation, CPU: Intel i7 @3.6 GHz, RAM: 8 GB).

Calibration Method
Step 1 Step 2 Total EM 0.3549s 0.0012s 0.3561s PM 0.0161s 0.0014s 0.0175s In summary, three conclusions are drawn: (1) using the proposed calibration method, all of the magnetometer errors are compensated, including the hard iron effect, soft iron effect, zero bias error, scale factor error, non-orthogonal error, and misalignment error; (2) using the standalone ellipsoid-fitting calibration algorithm, the total scale factor error and the total zero bias error are compensated, but the total misalignment error is uncompensated; and (3) the computational requirements of the proposed method is much lower than the computational requirements of the ellipsoid-fitting method.

Experiment
In this experiment, an SBG IG-500N MARG unit is used. As specified in the manual, the norm of the tri-axis magnetometer output vector is approximately equal to 1 with an arbitrary unit. The noise characteristics of the magnetometer and gyroscope are tested first. The magnetometer output noise is about 0.006 (1σ). The gyroscope output noise and bias are about 0.07 • /s (1σ) and 0.5 • /s, respectively, and the gyro bias errors are not compensated in advance. The sample frequency is set to be 100 Hz.
The unit is mounted to a cuboid frame as shown in Figure 9. The number of turns n is set to 2, and the interpolation angular interval ∆θ is set to 1 • . The frame is put on a flat wooden table, and rotated by hand over n cycles for each side of the frame. The simulation program is repeated ten times under both of the error conditions using Matlab, and the average computation time of these two algorithms is shown in Table 4. The step 1 in the EM requires large amounts of computation. The linear operations in the PM requires small computational amounts. The EM requires computations amounting to more than 10 times that of the PM. In summary, three conclusions are drawn: (1) using the proposed calibration method, all of the magnetometer errors are compensated, including the hard iron effect, soft iron effect, zero bias error, scale factor error, non-orthogonal error, and misalignment error; (2) using the standalone ellipsoidfitting calibration algorithm, the total scale factor error and the total zero bias error are compensated, but the total misalignment error is uncompensated; and (3) the computational requirements of the proposed method is much lower than the computational requirements of the ellipsoid-fitting method.

Experiment
In this experiment, an SBG IG-500N MARG unit is used. As specified in the manual, the norm of the tri-axis magnetometer output vector is approximately equal to 1 with an arbitrary unit. The noise characteristics of the magnetometer and gyroscope are tested first. The magnetometer output noise is about 0.006 (1σ). The gyroscope output noise and bias are about 0.07°/s (1σ) and 0.5°/s, respectively, and the gyro bias errors are not compensated in advance. The sample frequency is set to be 100 Hz.
The unit is mounted to a cuboid frame as shown in Figure 9. The number of turns n is set to 2, and the interpolation angular interval Δθ is set to 1°. The frame is put on a flat wooden table, and rotated by hand over n cycles for each side of the frame. By applying the PM, the total error coefficient matrix and total bias error vector are calculated as follows: By applying the EM, the total error coefficient matrix and total bias error vector are calculated as follows: The raw and calibrated tri-axis magnetometer output data during the calibration process are plotted in Figure 10. After calibration using the PM, the data lie on the circles parallel to the sides of the cuboid, centered at (0, 0). After calibration using the EM, the data lie on the circles unparalleled to the sides of cuboid, which means that the misalignment error still exists. The results of the experiment are consistent with the results of the simulation.
The raw and calibrated tri-axis magnetometer output data during the calibration process are plotted in Figure 10. After calibration using the PM, the data lie on the circles parallel to the sides of the cuboid, centered at (0, 0). After calibration using the EM, the data lie on the circles unparalleled to the sides of cuboid, which means that the misalignment error still exists. The results of the experiment are consistent with the results of the simulation. For further verification of the quality of the calibrations, the frame is installed on a turntable. By choosing the installation location, the magnetic disturbance from the turntable is reduced to a low level. The turntable is rotated over 1 cycle. Assuming the output angle from turntable as the ideal heading angle, the angle errors after calibration using the EM and PM are derived under the conditions of side 1 upwards, side 3 upwards and side 5 upwards, respectively, as shown in Figure 11. For further verification of the quality of the calibrations, the frame is installed on a turntable. By choosing the installation location, the magnetic disturbance from the turntable is reduced to a low level. The turntable is rotated over 1 cycle. Assuming the output angle from turntable as the ideal heading angle, the angle errors after calibration using the EM and PM are derived under the conditions of side 1 upwards, side 3 upwards and side 5 upwards, respectively, as shown in Figure 11. As shown in Table 5, by applying the EM, the angle errors are narrowed to 2.5° (1σ) with the maximum less than 4.2°. By applying the PM, the angle errors are narrowed to 0.5° (1σ) with the maximum less than 1.7°. The average computation time is shown in Table 6. The EM requires computations amounting to more than 10 times that of the PM, which is consistent with the simulation results.

Conclusions
This article proposed a complete tri-axis magnetometer calibration algorithm with a gyroscope auxiliary. A non-magnetic cuboid frame was required to mount the sensor set, and was rotated by hand during the calibration, which was inexpensive and easy to operate.
The calibration errors introduced by the gyro bias error, gyro misalignment error, cuboid frame non-orthogonality error, and rotation deviation error were analyzed, respectively. The results of the analyses indicated that the calibration errors were slight and acceptable. As shown in Table 5, by applying the EM, the angle errors are narrowed to 2.5 • (1σ) with the maximum less than 4.2 • . By applying the PM, the angle errors are narrowed to 0.5 • (1σ) with the maximum less than 1.7 • . The average computation time is shown in Table 6. The EM requires computations amounting to more than 10 times that of the PM, which is consistent with the simulation results. Table 6. Computation time (experiment, CPU: Intel i7 @3.6 GHz, RAM: 8 GB).

Conclusions
This article proposed a complete tri-axis magnetometer calibration algorithm with a gyroscope auxiliary. A non-magnetic cuboid frame was required to mount the sensor set, and was rotated by hand during the calibration, which was inexpensive and easy to operate. The calibration errors introduced by the gyro bias error, gyro misalignment error, cuboid frame non-orthogonality error, and rotation deviation error were analyzed, respectively. The results of the analyses indicated that the calibration errors were slight and acceptable.
The computation time of the proposed method was compared with the ellipsoid-fitting calibration method, which indicated that the linear operations provided this algorithm with low computational complexity.
The calibration results of the ellipsoid-fitting method and the proposed method were compared using simulations and experiments. The tri-axis magnetometer output errors were calibrated using the proposed method, including the hard iron effect, soft iron effect, zero bias error, scale factor error, non-orthogonal error, and misalignment error, while the total misalignment error, introduced by the misalignment error, soft iron error, and non-orthogonality error, was unable to be calibrated under the stand-alone ellipsoid fitting calibration algorithms. The results from the simulations and experiments were consistent. After calibration using the proposed method, the heading error was reduced to 0.5 • (1σ).