Next Article in Journal
Intuitive Cognition-Based Method for Generating Speech Using Hand Gestures
Next Article in Special Issue
Hybrid Magnetic–Inductive Angular Sensor with 360° Range and Stray-Field Immunity
Previous Article in Journal
Detection of Diabetic Eye Disease from Retinal Images Using a Deep Learning Based CenterNet Model
Previous Article in Special Issue
Low Field Optimization of a Non-Contacting High-Sensitivity GMR-Based DC/AC Current Sensor
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Review

Magnetic Field Sensors’ Calibration: Algorithms’ Overview and Comparison

by
Konstantinos Papafotis
*,
Dimitris Nikitas
and
Paul P. Sotiriadis
*
Department of Electrical and Computer Engineering, National Technical University of Athens, 157 80 Athens, Greece
*
Authors to whom correspondence should be addressed.
Sensors 2021, 21(16), 5288; https://doi.org/10.3390/s21165288
Submission received: 7 June 2021 / Revised: 17 July 2021 / Accepted: 22 July 2021 / Published: 5 August 2021
(This article belongs to the Special Issue Magnetic Sensors 2021)

Abstract

:
The calibration of three-axis magnetic field sensors is reviewed. Seven representative algorithms for in-situ calibration of magnetic field sensors without requiring any special piece of equipment are reviewed. The algorithms are presented in a user friendly, directly applicable step-by-step form, and are compared in terms of accuracy, computational efficiency and robustness using both real sensors’ data and artificial data with known sensor’s measurement distortion.

1. Introduction

Magnetic field sensors (magnetometers) are nowadays widely used in a plethora of commercial, industrial, marine, aerospace and military applications. Their applications include but not limited to navigation and attitude estimation, geophysical surveys, archaeology, entertainment devices, consumer electronics and others.
In most applications, sensor’s calibration is essential in order to achieve the desirable accuracy level. The purpose of magnetic field sensors’ calibration is a twofold. First, as in the case of every measurement unit, calibration ensures that the measurement of the standalone sensor corresponds to the actual value of the magnetic field. To do so, calibration must compensate for all static (manufacturing imperfections etc.) and active (temperature, humidity, etc.) phenomena effecting the accuracy of the sensor’s measurement. In addition, when a magnetic sensor is embedded in a larger system, other components of the system may cause disturbances (both static and active ones) to the local magnetic field. Static disturbances are usually caused by magnetic and ferromagnetic materials in the vicinity of the sensor; called hard-iron distortion and soft-iron distortion respectively (more information are given in Section 2). Mechanical or electronic structures embedded in the system, such as motors and coils could also actively distort the local magnetic field and cause significant measurement error.
This review paper focuses on algorithms correcting the dominant linear time-invariant (static) measurement errors, requiring no special piece of equipment for their application. Such algorithms are most commonly used for in-situ calibration of magnetic field sensors which are usually in chip form and embedded in larger systems. The paper presents seven representative calibration algorithms for three-axis magnetometers and compares them in terms of accuracy, robustness, computational efficiency and ease of deployment. The seven algorithms are briefly presented, to introduce all required mathematical expressions, and are summarized in an easy-to-develop, step-by-step form. For the details of the algorithms, the reader is referred to the original works.
The selection of the particular algorithms was done based on their popularity and on our attempt to present as many different calibration approaches as possible. The TWOSTEP [1] algorithm is one of the first algorithms that addressed the full calibration problem (and probably the most popular one). At a later time, Elkaim and Vasconcelos [2] proposed a geometric approach of TWOSTEP which is also very popular. At the same time, Dorveaux et al. [3] offered a non-linear formulation of the problem and they treated it in an innovative, strictly iterative way. In addition, Wu and Shi [4] suggested the most complete formulation of the calibration problem as an optimal maximum likelihood estimation one. The TWOSTEP algorithm, as well as the algorithms proposed by Vasconcelos et al. and Wu et al., consist of a first step deriving an initial solution, and, a second step for improving it. On the other hand, Papafotis and Sotiriadis [5] recommended an iterative approach based on a twofold minimization, which was shown to be extremely effective. Furthermore, a real-time approach by Crassidis et al. [6] using the popular Kalman Filter is discussed. Finally, to represent the recent trends towards Machine Learning, an Artificial Intelligence (AI) method applying Particle Swarm Optimization on the estimation problem is explored [7].
Please note that this review focuses on works for in-situ calibration of three-axis magnetic field sensor without using any special piece of equipment or any other additional sensor. Thus, several interesting works dealing with magnetometer’s calibration, in combination with inertial sensors, [8,9,10,11,12] are not included in this work.
The rest of the paper is organized as follows. First, a standard error model for three-axis magnetic field sensors is presented in Section 2. In Section 3, Section 4, Section 5, Section 6, Section 7, Section 8 and Section 9, seven representative algorithms are discussed in chronological order of publication. In Section 10, a method for generating artificial data is proposed and algorithms are evaluated via extensive Monte Carlo simulation to identify their performance. In addition, the algorithms are evaluated using several real sensor’s measurements in order to evaluate their performance under real-world conditions. Finally, Section 11 summarizes our findings and provides brief comments for each algorithm. The notation used along the paper is presented in Table 1.

2. Magnetic Field Sensor’s Error Sources and Measurement Model

In this section, the most important linear, time-invariant error sources of three-axis magnetic field sensors are presented. Based on them, a mathematical model relating the sensor’s measurement with the actual value of the magnetic field is derived.
The total output error of a magnetic sensor is a combination of several error sources related to the sensing element itself, the instrumentation electronics, manufacturing imperfections and distortions caused by magnetic and ferromagnetic materials in the vicinity of the sensor. The linear, time-invariant error sources with the most significant contribution in the total sensor’s error, are listed below:
  • Bias, or offset; all magnetic sensors suffer from bias, which is a constant distortion. In many cases, it is the most important defect in the sensor’s overall error. A  3 × 1 vector, h s , is used to model it.
  • Scale-factor error represents the input-output gain error of the sensor. It is modeled by a 3 × 3 diagonal matrix, T s f .
  • Cross-coupling or non-orthogonality inaccuracies are resulted by the non-ideal alignment of the sensor’s axes during manufacturing and are modeled by a 3 × 3 matrix, T c c .
  • Soft-iron distortion is caused by ferromagnetic materials in the vicinity of the sensor, attached to the sensor’s coordinate frame. Those materials do not generate their own magnetic field, but instead alter the existing magnetic field locally, resulting in a measurement discrepancy. This effect is modeled by a 3 × 3 matrix, T s i .
  • Hard-iron distortion is due to magnetic materials attached to the sensor’s coordinate frame. As a consequence of the persistent magnetic field created by those materials, the sensor’s output has a constant bias. Hard-iron distortion is modeled by a 3 × 1 vector, h h i .
  • Random noise is the stochastic error in the sensor’s output. It is induced by the sensor’s mechanical and electrical architecture. It is modeled by a 3 × 1 vector, ε , and it is most commonly assumed to be a sequence of white noise, i.e.,  ε N ( 0 , σ 2 ) .
Let m be the 3 × 1 true magnetic field vector and y be the 3 × 1 measurement vector. With the aforementioned error terms in mind, a widely accepted and well-referenced measurement model for a three-axis magnetometer is the following [1,2,4,5,6,7,13,14]
y = T s f T c c T s i m + h h i + h s + ε
In most applications, the exact contribution of each error term in (1) is of no concern and, thus, instead of (1), most calibration algorithms use the following, compact form of (1)
y = T m + h + ε
where T T s f T c c T s i and h T s f T c c h h i + h s .
This work focuses on algorithms intended to be used with magnetic field sensors requiring no special piece of equipment. In such cases, the calibration is done in the sensor’s (body) coordinate frame implying that both the measurement vector, y and the true magnetic field vector, m in (2) are expressed in the senor’s coordinate frame.
Note that when expensive laboratory equipment is not available, both the calibration parameters T and h in (2), and the magnetic field vector, m, are unknown. Thus, in most works, multiple measurements of the local (Earth’s) magnetic field are used to derive T and h. Note that the Earth’s magnetic field varies with location and time and its value (magnitude and direction) is only approximately known by magnetic models such as International Geomagnetic Reference Field model (IGRF) [15]. However it is reasonable to assume that the magnitude of the magnetic field is (locally) constant during the calibration procedure. Based on this fact, most authors formulate an optimization or an estimation problem to derive T and h.

3. Alonso and Shuster (TWOSTEP)

The TWOSTEP algorithm [1] consists of an analytical centering approach [16,17] for its first step, while in the second step the solution is optimized numerically. The authors initially solved the problem of bias, h, determination when attitude is not known [18] and then extended their method to determine matrix T as well [1].
TWOSTEP is motivated by the assumption that matrix T should not be far from a pure rotation. Therefore by applying polar decomposition it can be written as T = ( I 3 × 3 + D ) 1 O where O is an orthogonal matrix and D is a symmetric 3 × 3 matrix so as ( I 3 × 3 + D ) 1 to be positive definite. Matrix O can be integrated into vector m since it does not alter its norm. The equivalent measurement model is
y = T ^ m ^ + h + ε
where
T ^ ( I 3 × 3 + D ) 1 m ^ O m .
Therefore, for the full calibration, D and h must be estimated. To this purpose, a set of N measurements, y k , k = 1 , 2 , N , is used and the corresponding effective measurements z k , k = 1 , 2 , N , are defined as
z k y k 2 m k ^ 2 = y k 2 m k 2 .
The last ones can be decomposed into a deterministic part plus an approximately Gaussian noise term, υ k with mean μ k and variance σ k 2 , i.e., υ k N ( μ k , σ k 2 ) , given by
μ k = 3 σ 2
σ k 2 = 4 σ 2 ( ( I 3 × 3 + D ) y k h ) T ( ( I 3 × 3 + D ) y k h ) + 6 σ 4 .
Since D and h are unknown, the variance σ k 2 is assumed to be similar to measurement’s output error variance σ 2 . Hence μ k and σ k 2 can be assumed independent of k. To estimate D and h, Alonso and Shuster define the auxiliary quantities E D 2 + 2 D and c ( I + D ) h , as well as the estimation vector θ containing the elements of the 3 × 1 vector c and those of the 3 × 3 symmetric matrix E, which is formed as θ = [ c T E 11 E 22 E 33 E 12 E 13 E 23 ] T .
TWOSTEP algorithm functions on the estimation vector θ and thus on the auxiliary parameters, E and c and not on the actual calibration parameters, D and h. The transformation from E and c back to D and h is described in (13) and (14).

3.1. Initial Estimate

For every measurement, y k , k = 1 , 2 , N , the following auxiliary variables are defined
S k = [ y k , 1 2 y k , 2 2 y k , 3 2 2 y k , 1 y k , 2 2 y k , 1 y k , 3 2 y k , 2 y k , 3 ]
L k = [ 2 y k T S k ] .
The centering approximation is done using the following weighted averages along with their corresponding centered values
z ¯ σ ¯ 2 k = 1 N 1 σ k 2 z k z ˜ k z k z ¯
L ¯ σ ¯ 2 k = 1 N 1 σ k 2 L k L ˜ k L k L ¯
μ ¯ σ ¯ 2 k = 1 N 1 σ k 2 μ k μ ˜ k μ k μ ¯
where
σ ¯ 2 k = 1 N 1 σ k 2 1 .
The first estimation of θ is the centered one given by
θ ˜ = P ˜ θ θ k = 1 N 1 σ k 2 z ˜ k μ ˜ k L ˜ k T
P ˜ θ θ 1 = k = 1 N 1 σ k 2 L ˜ k T L ˜ k
with P ˜ θ θ denoting the centered covariance matrix and F ˜ θ θ denoting the centered Fischer information matrix.

3.2. Solution Improvement Step

The second step improves the previous estimate of vector θ , derived in (8), via Gauss-Newton method using the centered estimate θ ˜ as the initial guess. The estimation is updated as follows
θ i + 1 = θ i F ˜ θ θ + 1 σ ¯ 2 L ¯ ϕ ( θ i ) T L ¯ ϕ ( θ i ) 1 g ( θ i )
where
v = ( I 3 × 3 + E ) 1 c
ϕ ( θ ) = [ 2 v T v 1 2 v 2 2 v 3 2 2 v 1 v 2 2 v 1 v 3 2 v 2 v 3 ]
g ( θ ) = P ˜ θ θ 1 ( θ θ ˜ ) 1 σ ¯ 2 z ¯ L ¯ θ + c T v μ ¯ L ¯ T ϕ ( θ )
with v j denoting the j t h element of vector v. At every iteration the 3 × 3 symmetric matrix E and the 3 × 1 vector c are updated according to the current estimation vector θ i using
c = θ 1 θ 2 θ 3 and E = θ 4 θ 7 θ 8 θ 7 θ 5 θ 9 θ 8 θ 9 θ 6 .
Alonso and Shuster define the following quantity in order to establish a stop condition for the Gauss–Newton method.
η i θ i + 1 θ i T F ˜ θ θ + 1 σ ¯ 2 L ¯ ϕ ( θ i ) T L ¯ ϕ ( θ i ) θ i + 1 θ i .
The iterations continue until η i became smaller than a predetermined threshold.
After sufficiently many iterations, an optimal estimation of matrix E * and of vector c * is derived. The derived solution is then used to find D * and h * . To this end we apply SVD [19] to the symmetric matrix E * , i.e.,
E * = U S U T
where S = diag ( s 1 , s 2 , s 3 ) , U O ( 3 ) . Then find the diagonal matrix W = diag ( w 1 , w 2 , w 3 ) satisfying S = 2 W + W 2 . Typically, the elements of S are much smaller than unity [1] so a real solution exists with the diagonal entries of W being w j = 1 + 1 + s j , j = 1 , 2 , 3 .
Combining the above, the estimates of matrix D * and bias vector h * are given by
D * = U W U T
h * = ( I 3 × 3 + D * ) 1 c *
and are related to the calibration parameters T and h of the measurement model (2) as follows
T = ( I 3 × 3 + D * ) 1 and h = h * .
Summarizing, when the centered estimation is near the ground truth value the Gauss–Newton method typically converges rapidly. The authors verified the robustness of their method via simulations assuming either white noise or colored noise. TWOSTEP is also suitable for on-orbit calibration using IGRF data [15]. The algorithm is summarized in Algorithm 1.
Algorithm 1: Alonso and Shuster (TWOSTEP) [1]
Step 1: Calculate z k , L k , for k = 1 , 2 , , N
        by using (4)–(6)  
Step 2: Calculate the centered values z ˜ k , L ˜ k for k = 1 , 2 , , N (7)
Step 3: Calculate centered estimate θ ˜ and covariance matrix P ˜ θ θ (8)
Step 4: Extract c and E from θ following (11)
Step 5: Calculate ϕ ( θ ) and g ( θ ) from (10)
Step 6: Update θ using (9)
Step 7: Calculate η following (12)
Step 8: Repeat steps 4–7 until η is sufficiently small  
Step 9: Apply SVD on E * (13) and define matrix W
Step 10: Calculate D * , h * (14) and T , h (15)

4. Crassidis et al.

The authors of [6] were motivated by the fact that real-time applications demand real-time calibration methods. To this end, based on the problem formulation (3) established in TWOSTEP [1], Crassidis et al. formulate a real-time identification problem for the derivation of the calibration parameters D and h and solve it using the extended Kalman Filter (EKF) approach. Note that the authors of [6] have proposed two more algorithms dealing with the online calibration of a magnetic field sensor in [20,21]. However, in this work we focus on the EKF based one presented in [6] which is the most efficient and popular one.
Following the problem formulation of TWOSTEP, the bias vector h and the symmetric matrix D are desired. The estimation vector θ is defined differently and contains h and D, structured as follows
θ = [ h T D 11 D 22 D 33 D 12 D 13 D 23 ] T .
Because the vector θ is constant, the state model is given by θ ˙ = 0 . The effective measurement is given by z k = y k 2 m k 2 (4) while the measurement’s model is given by z k = ϕ ( θ k ) + υ k where
ϕ ( θ k ) = y k T 2 D k + D k 2 y k + 2 y k T I 3 × 3 + D k h k h k 2
and effective measurement’s noise υ k N ( μ k , σ k 2 ) follows (5). At each iteration D k and h k are extracted from θ k according to (16). The propagation is as it follows
θ k + 1 = θ k + K k z k ϕ ( θ k )
P k + 1 = I 9 × 9 K k H ( θ k ) P k
K k = P k H T ( θ k ) H ( θ k ) P k H T ( θ k ) + σ k 2 1
where P k is the covariance of the estimated parameters for h and D at step k. The matrix H ( θ k ) is the linearization matrix of ϕ ( θ k ) and is defined as  
H ( θ k ) = [ 2 y k T ( I 3 × 3 + D k ) 2 h k T S k F k + 2 J k ]
where
S k = [ y k , 1 2 y k , 2 2 y k , 3 2 2 y k , 1 y k , 2 2 y k , 1 y k , 3 2 y k , 2 y k , 3 ]
J k = [ y k , 1 h k , 1 y k , 2 h k , 2 y k , 3 h k , 3 y k , 1 h k , 2 + y k , 2 h k , 1 y k , 1 h k , 3 + y k , 3 h k , 1 y k , 2 h k , 3 + y k , 3 h k , 2 ]
F k = Δ 1 0 0 2 D k , 12 2 D k , 13 0 0 Δ 2 0 2 D k , 12 0 2 D k , 23 0 0 Δ 3 0 2 D k , 13 2 D k , 23 D k , 12 D k , 12 0 Δ 4 D k , 23 D k , 13 D k , 13 0 D k , 13 D k , 23 Δ 5 D k , 12 0 D k , 23 D k , 23 D k , 13 D k , 12 Δ 6
and
Δ 1 = 2 ( 1 + D k , 11 ) Δ 2 = 2 ( 1 + D k , 22 ) Δ 3 = 2 ( 1 + D k , 33 ) Δ 4 = 2 + D k , 11 + D k , 22 Δ 5 = 2 + D k , 11 + D k , 33 Δ 6 = 2 + D k , 22 + D k , 33 .
  The noise variance of the measurements, σ k 2 , can be assumed to be constant and equal to σ 2 as in TWOSTEP. Given a set of N measurements, the EKF provides an optimal estimation vector θ * = θ N from which an optimal vector h * = h N and a matrix D * = D N can be extracted according to (16). Therefore, the calibration parameters (2) are given by
T = ( I 3 × 3 + D * ) 1 and h = h * .
Under certain conditions, e.g., fast changing data, this approach of sequential calibration may have some advantages in terms of computational complexity and adaptation. The algorithm is summarized in Algorithm 2.
Algorithm 2: Crassidis et al. [6]
Step 1: Initialize θ and k = 0
Step 2: for each measurement do:
                Calculate z k (4)
                Extract D k and h k from θ k (16)
                Calculate S k , J k , F k (20) and H ( θ k ) (19)
                Calculate Kalman Gain K k (18)
                Update estimation: θ k θ k + 1
                Update covariance matrix: P k P k + 1 (18)
                k k + 1
Step 3: Extract D * and h * from θ * (16)
Step 4: Calculate T and h (22)

5. Dorveaux et al.

An iterative algorithm for the calibration of magnetic field sensors based on iterations of a least-squares problem is introduced in [3]. In the beginning of the algorithm, the measurements lie on an ellipsoid according to (2). In each iteration, the measurements move from the initial ellipsoid to the unit sphere, following a cost function minimization algorithm.
The authors in [3] use the following variation of the measurement model of (2)
m = A y + B
where A = T 1 , B = T 1 h and the measurement noise, ε , is neglected.
The algorithm begins by considering an initial estimate of the magnetic field vectors, denoted by m ˜ k ( 0 ) and defined as
m ˜ k ( 0 ) = y k , k = 1 , 2 , , K .
In every iteration, the following cost function is formulated and minimized using the least squares method.
J ( A , B , n ) = k = 1 K A m ˜ k ( n ) + B m ˜ k ( n ) m ˜ k ( n ) 2
where n = 1 , 2 , , N denotes the n t h iteration. Let A n and B n be the resulting matrices from the minimization of (25). Every iteration ends with using A n and B n to update the estimates of the magnetic field vectors as
m ˜ k ( n + 1 ) = A n m ˜ k ( n ) + B n , k = 1 , 2 , , K .
From (26) we can express the magnetic field estimates m ˜ k ( n ) using the measurement vectors y k as
m ˜ k ( n ) = A ˜ n y k + B ˜ n , k = 1 , 2 , , K
where A ˜ n and B ˜ n are iteratively defined as
A ˜ n = A n A ˜ n 1 and B ˜ n = A n B ˜ n 1 + B n .
To determine when the algorithm has reached an acceptable solution, we define the following cost
J s t o p ( A n , B n ) = B n + A n I 3 × 3 .
The iterations stop when J s t o p is sufficiently small. Note that the original manuscript does not provide an explicit condition to stop iterations. However it is reasonable to terminate the algorithm when contribution of the updated A n and B n to the calibration parameters A ˜ n and B ˜ n is negligible (see (28)). The estimate m ˜ ( N ) derived at the Nth iteration represents the calibrated data and it is:
m k = m ˜ k ( N ) , k = 1 , 2 , , K
The derived matrices A ˜ N and B ˜ N are related to the calibration parameters T and h of the measurement model (2) as follows
T = A ˜ N 1 and h = A ˜ N 1 B ˜ N .
Finally, the estimates m ˜ k ( N ) , k = 1 , 2 , , K , derived at the N t h iteration represent the calibrated measurement vectors. The algorithm is summarized in Algorithm 3.
Algorithm 3: Dorveaux et al. [3]
Step 1: Initialize m ˜ k ( 0 ) using (24).
Step 2: Minimize (25) using least squares and derive A n and B n .
Step 3: Use A n and B n to calculate m ˜ k ( n + 1 ) from (26).
Step 4: Calculate A ˜ n and B ˜ n using (28).
Step 5: Evaluate the cost function J s t o p ( A n , B n ) from (29).
Step 6: Repeat steps 2-5 until J s t o p is sufficiently small.
Step 7: Use A ˜ N and B ˜ N to calculate T and h using (31).

6. Vasconcelos et al.

The authors of [2] consider that magnetometers’ measurements lie on a ellipsoid manifold following the measurement model (2). First, they derive an initial estimate of the calibration parameters T and h by finding the ellipsoid that fits best to the given data. Then, they use the measurement model of (2) to formulate a maximum likelihood estimation problem and derive an improved estimate of the calibration parameters T and h.
From (2), the magnetic field vector is expressed as m = T 1 ( y h ) T 1 ε . Assuming that the magnitude of the magnetic field is constant during the calibration procedure we can write the following unconstrained optimization problem to derive T and h
minimize T , h k = 1 K T 1 ( y k h ) 1 σ k 2 .
Here σ k denotes the standard deviation of the measurement noise in the k t h measurement, assuming it is the same for all three axes and equal to σ . Without loss of generality, the magnitude of the magnetic field is assumed to be equal to one. A possible relaxation of this soft assumption is provided by Springmann [22] who addresses the problem of time-varying bias. To solve (32), the authors define the following cost function and then minimize it using the Newton’s method
J ( x ) k = 1 K T ^ ( y k h ) 1 σ k 2
where T ^ = T 1 and
x = vec ( T ^ ) T h T T .
The vector x is updated in every Newton’s iteration as follows
x ( + ) = x ( ) 2 J ( x ) | x = x ( ) 1 J ( x ) | x = x ( )
where J ( x ) is the gradient vector and 2 J ( x ) is the Hessian matrix of the cost function. For both J ( x ) and 2 J ( x ) , the authors in [2] provide analytical expressions which are presented in Appendix A.1.

Initial Estimate

Solving (32) using the Newton’s method requires a good initial estimate of the calibration parameters, T ^ and h. Vasconcelos et al. use a previous work on nonlinear estimators for strapdown magnetometers by Foster and Elkaim [23,24], to derive a good initial estimate. Solving the ellipsoid equation m k = T 1 ( y k h ) = 1 for every k is equivalent to solving the following pseudo-linear least squares estimation problem by re-arranging the terms as follows
L p = b
where, by writing each measurement vector as y k = y k x y k y y k z T , k = 1 , 2 , , K , it is
L = y 1 x 2 y 1 x y 1 y y 1 x y 1 z y 1 y 2 y 1 y y 1 z y 1 x y 1 y y 1 z 1 y K x 2 y K x y K y y K x y K z y K y 2 y K y y K z y K x y K y y K z 1
and
b = y 1 z 2 y 2 z 2 y K z 2 T .
The vector p is derived as
p = A B C D E G H I J T = ( L T L ) 1 L T b .
The initial estimates of the calibration parameters are derived as
T ^ ( 0 ) = 1 α 0 0 1 α tan ( ρ ) 1 b sec ( ρ ) 0 1 α tan ( ρ ) tan ( λ ) sec ( ϕ ) tan ( ϕ ) 1 b sec ( ρ ) tan ( λ ) sec ( ϕ ) 1 c sec ( λ ) sec ( ϕ )
and
h ( 0 ) = 1 2 α 1 β 1 β 2 β 3 T
where
a = 1 2 α 1 ( 4 D + E 2 ) α 2 1 / 2 b = 1 2 α 1 ( 4 A + C 2 ) α 2 1 / 2 c = 1 2 α 1 ( 4 D A B 2 ) α 2 1 / 2 tan ( ρ ) = 1 2 α 1 ( 2 B + E C ) ( α 1 ) 1 / 2 tan ( ϕ ) = ( B E 2 C D ) ( α 1 ) 1 / 2 tan ( λ ) = E ( α 1 α 3 1 ) 1 / 2
and
β 1 = 2 B H + B E I 2 C D I 4 D G + E C H E 2 G β 2 = 2 A E I + 4 A H B C I 2 B G + C 2 H C E G β 3 = 4 D I A 2 D G C + E G B I B 2 2 E H A + C B H .
The auxiliary variables α 1 , α 2 and α 3 are defined as
α 1 = B 2 + D C 2 + 4 D A + A E 2 B E C α 2 = 4 A E 2 J E 2 G 2 4 B E C J + 2 E C H G + 2 B E I G 4 E H A I 4 D I C G C 2 H 2 + 4 D A I 2 + 2 C B H I 4 D G 2 + 4 D C 2 J + 4 B H G 4 A H 2 B 2 I 2 4 B 2 J + 16 D A J α 3 = E 4 A C B E 3 + E 2 C 2 D 2 B 2 E 2 + 8 D A E 2 4 D B 2 + 16 D 2 A .
One contribution of Vasconcelos et al., advancing the existing initial step approach suggested in [23], was the derivation of the aforementioned explicit and non-trivial expressions. In addition, Vasconcelos et al. state that their proposed algorithm is applicable even when the magnitude of the magnetic field is not constant during the measurement, similarly to TWOSTEP and Crassidis et al. algorithm [6]. The algorithm is summarized in Algorithm 4.
Algorithm 4: Vasconcelos et al. [2]
Initial Estimate
Step 1: Use the sensors’ measurements y k , k = 1 , 2 , , K and form A and b according to (37) and (38), respectively.
Step 2: Calculate p using (39)
Step 3: Derive the initial estimates T ^ ( 0 ) and h ( 0 ) using (40) and (41), respectively.
Newton Method
Step 4: Use the initial estimates T ^ ( 0 ) and h ( 0 ) to initialize x according to (34).
Step 5: Update x using (35).
Step 6: Evaluate the cost function J ( x ) of (33).
Step 7: Repeat Steps 5–6 until J ( x ) becomes sufficiently small.
Step 8: Split x into T ^ and h and calculate T = T ^ 1 .

7. Ali et al.

The authors in [7] propose a Particle Swarm Optimization (PSO) [25] - based calibration algorithm that estimates the bias, the scale and nonorthogonality factors. The main advantage of this algorithm is its simplicity of implementation since the optimization is heuristic and does not depend on calculation of gradients, unlike other optimization techniques mentioned in this paper. It can be classified as an AI [26] approach.
The authors in [7] use (2) and a set of N sensor’s measurements to form the following optimization problem for deriving the calibration parameters T and h
minimize T , h k = 0 N y k 2 m k 2 2 .
where J k = 0 N y k 2 m k 2 2 is called the fitness value.
Function J depends on T and h which are conveniently combined into the single vector x R 12 ,
x = h vec ( T T ) .
For a swarm of S particles, the position x i R 12 and the velocity v i R 12 of the i t h particle can be computed using [25]
v i k = v i k 1 + c 1 r 1 i k 1 p i k 1 x i k 1 + c 2 r 2 i k 1 p g k 1 x i k 1
x i k = x i k 1 + v i k
for i = 1 , 2 , , S where k denotes the new value while k 1 the old value. Also p i denotes the i t h ’s particle best position, p g denotes the swarm’s best position, c 1 and c 2 are the acceleration coefficients, w is the inertial weight factor and r 1 i , r 2 i are random numbers uniformly distributed within the range [ 0 , 1 ] . Typical values of these quantities are c 1 = c 2 = 2 , w = 1 and the number of particles S is usually between 20 and 65.
Therefore, at each iteration k, each particle’s fitness value J ( x i k ) is calculated and quantities p i and p g are updated accordingly. The authors suggest three different stop criteria. Specifically, the iterations stop either when the fitness value J of a particle is smaller than a predetermined threshold, or after a maximum number of iterations, or when the change of J becomes insignificant with iterations. Upon termination of the algorithm, parameters T and h (2) are extracted from the swarms’s optimal solution p g according to
h vec ( T T ) = p g .
Following the general concept of applying AI optimization algorithms, as was introduced in [7], one can also consider using more modern versions of the standard PSO, e.g., [27,28,29]. They are typically found as built-in functions in computational suites such as MATLAB [30]. The algorithm is summarized in Algorithm 5.
Algorithm 5: Ali et al. [7]
Step 1: Initialize x i , v i for i = 1 , 2 , , S
        and set p i = x i
Step 2: Find j = i | i = 1 , 2 , , S and J ( p i ) min
        Particle i best: J m i n i J ( p i )
        Global best: p g p j and J m i n J ( p j )
Step 3: for each particle i do
                Update x i , v i (47)
                Calculate J ( x i ) (45)
                if  J ( x i ) < J m i n i
                   J m i n i J ( x i ) and p i x i
                   if  J ( x i ) < J m i n
                       J m i n J ( x i ) and p g x i
Step 4: Repeat Step 3 until an exit condition is met
Step 5: Extract T and h from p g (48)

8. Wu and Shi

The authors of [4], formulate the calibration of a three-axis magnetometer as a maximum likelihood estimation problem which is solved using the Gauss-Newton method.
Starting from the measurement model of (2), Wu and Shi observed that by considering the QR decomposition T 1 = Q R , where Q O ( 3 ) and R U ( 3 ) , (2) is written as
y = R 1 Q T m + h + ε .
Defining m ^ Q T m , we observe that m ^ = m since Q O ( 3 ) . Also setting T ^ R 1 we have that
y = T ^ m ^ + h + ε .
Using the above transformation, the authors reduce the unknown model parameter variables from 12 (9 for T and 3 for h) to 9 (6 for R since R is upper triangular and 3 for h). Note that using (50), the calibration procedure now aims at finding the calibration parameters T ^ and h while the magnetic field vector m ^ is also unknown.
Using a set of K measurements and (50), the authors formulate the following maximum likelihood estimation problem
minimize T ^ , h , m ^ k k = 1 K y k T ^ m ^ k h 2 subject to m ^ k = 1 , k = 1 , 2 , , K .
Without loss of generality, the authors, constrained the magnitude of the magnetic field to be equal to one. Based on (51), the following Lagrange function is formulated
J ( x ) = k = 1 K y k T ^ m ^ k h 2 + λ k m ^ k 2 1
where
x = vec ( T ^ ) T , h T , m ^ 1 T , m ^ 2 T , , m ^ K T , λ 1 , λ 2 , , λ K T
and λ k , k = 1 , 2 , , K are positive Lagrange coefficients for the unit norm constrain. Note that since T ^ is an upper triangular matrix, the lower triangular elements of T ^ are excluded from x. The minimization of (52) and the estimation of x are done using the Gauss-Newton method as follows
x ( + ) = x ( ) 2 J ( x ) | x = x ( ) 1 J ( x ) | x = x ( )
where J ( x ) is the Jacobian vector and 2 J ( x ) is the Hessian matrix of the Lagrange function. For both J ( x ) and 2 J ( x ) , the authors provide analytical expressions which are presented in Appendix A.2.

Initial Estimate

Solving (51) using the Gauss-Newton method requires a good initial estimate of the unknowns. To find one, the authors of [4] use the unit magnitude constrain and the equation 1 = R y k h 2 which after some manipulation, is written as
y k T y k T y k T 1 vec ( A ) b c Y k z = 0 , k = 1 , 2 , , K
where A = R T R , b = 2 R T R h and c = h T R T R h . Defining Y = Y 1 T Y 2 T Y K T T , from (55) it is
Y z = 0
The authors, solve (56) using the least squares method and denote the solution z e = vec ( A e ) T b e T c e T = min Y z 2 . They derive z e as the eigenvector of Y T Y corresponding to its minimum (or zero) eigenvalue. Using z e , the vector z is derived as z = α z e , where α = 4 / b e T A e 1 b e 4 c e . Extracting vec ( A ) , b and c from z, the initial estimates of the unknowns, T ^ ( 0 ) , h ( 0 ) , m ^ k ( 0 ) and λ k ( 0 ) are defined as follows:
T ^ ( 0 ) = R 1 = chol ( A ) h ( 0 ) = A 1 b / 2 m ^ k ( 0 ) = T ^ ( 0 ) 1 ( y k h ) , k = 1 , 2 , , K λ k ( 0 ) = 0 , k = 1 , 2 , , K
where chol ( · ) is the Cholesky factorization.
Finally, an alternative version of Wu’s and Shi’s algorithm is proposed by Cao et al. in [13], where a different method for the initial estimate is presented, and the second step is identical. The algorithm is summarized in Algorithm 6.
Algorithm 6: Wu and Shi [4]
Initial Estimate
Step 1: Calculate Y k , k = 1 , 2 , , K from (55) and form the matrix Y = Y 1 T Y 2 T Y K T T .
Step 2: Find the eigenvector of Y T Y corresponding to its minimum (or zero) eigenvalue and denote it as z e = vec ( A e ) T b e T c e T .
Step 3: Calculate z = a z e where α = 4 / b e T A e 1 b e 4 c e .
Step 4: Extract vec ( A ) , b and c from z.
Step 5: Calculate an initial estimate of the unknowns using (57).
  Gauss–Newton Method
Step 6: Use the initial estimates to initialize the vector x of (53)
Step 7: Update x using (54).
Step 8: Evaluate the cost J ( x ) of (52).
Step 9: Repeat steps 7-8 until J ( x ) becomes sufficiently small.

9. Papafotis and Sotiriadis (MAG.I.C.AL.)

The authors in [5] use (2) and a set of K sensor’s measurements to form the following optimization problem for deriving the calibration parameters T and h
minimize T , h , m k k = 1 K y k T m k h 2 subject to m k = 1 , k = 1 , 2 , , K
where, without loss of generality, the magnitude of the magnetic field is constrained to be equal to one. In order to solve (58) they propose an iterative algorithm, based on the solution of a linear least-squares problem.
The algorithm begins by initializing the magnetic field vectors, m k , as 
m k = y k y k , k = 1 , 2 , , K
and rewriting (2) in a matrix form as follows:
Y = L G + E
where
Y = y 1 y 2 y K
L = T h
G = m 1 m 2 m K 1 1 1
E = ε 1 ε 2 ε K .
In every iteration, (60) is solved for L using the least squares method, minimizing the total squared error E T E F 2
L = Y G T ( G G T ) 1
From the calculated L, an updated set of calibration parameters T and h is extracted from (61b). Using them, the magnetic field vector is updated as
m k = m k ˜ m k ˜ , k = 1 , 2 , , K
where
m ˜ k = T 1 ( y k h ) , k = 1 , 2 , , K .
Every iteration ends by updating the matrix G using the updated vectors m k , k = 1 , 2 , , K . Iterations stop when a small value of the following cost function is achieved
J ( T , h ) = k = 1 K m ˜ k 2 1 2 .
MAG.I.C.AL. algorithm is summarized in Algorithm 7.   
Algorithm 7: Papafotis and Sotiriadis (MAG.I.C.AL.) [5]
Step 1: Initialize m k using (59).
Step 2: Calculate L using (62).
Step 3: Extract T and h from L using (61b).
Step 4: Update m k using (63) and (64) and use it to update G.
Step 5: Evaluate the cost-plus-penalty function J from (65).
Step 6: Repeat steps 2-5 until J ( T , h ) is sufficiently small.

10. Algorithm Evaluation and Comparison

In this Section, the performance of the presented algorithms are evaluated in terms of accuracy, robustness, and execution speed. Firstly, we evaluate the performance of the seven algorithms using multiple sets of synthetic data where the calibration parameters T and h, as well as the measurement noise characteristics are predefined and known. By doing so, we are able to accurately determine the algorithms’ accuracy and robustness. Then multiple datasets of two different low-cost magnetic field sensors are used to verify the algorithms’ performance under real-world conditions.

10.1. Synthetic Data Generation

We designed a procedure to generate synthetic data effectively, in order to examine each of the aforementioned algorithm’s performance across a range of noise variance and measurement sample size. The authors of TWOSTEP [18] propose a typical scenario of assuming the magnetic vector spinning with a constant angular velocity. On the other hand, Wu and Shi [4] suggest a specific sequence of 3D rotations using Euler Angles, applied on a constant known magnetic vector m. In the same page, Papafotis and Sotiriadis [5] recommend a sequence of 12 approximate orientations. Another alternative is to make use of a set of random, yet normalized, vector fields, which however demands a reasonable amount of samples.
Because none of the described algorithms guarantees that it will function properly under an arbitrary dataset, we propose an efficient method to span S O ( 3 ) , following [31], so as to provide the algorithms with substantial, non-redundant information and to compare them fairly. After extensive simulation, it was observed that the recommended method was very effective in spanning the 3D rotation space.
Our method’s effectiveness lies in distributing the points on the sphere m = 1 , more evenly by using the canonical Fibonacci lattice mapping [31,32]. Generating a Fibonacci sphere is an extremely fast and effective approximate method to evenly distribute points on a sphere.
This way S O ( 3 ) is sufficiently represented even with only a small dataset. An algorithm for generating K vectors distributed on a Fibonacci sphere is presented in detail in Algorithm 8.
Considering K vectors, m k , k = 1 , 2 , , K distributed on a Fibonacci sphere, we continue with generating matrix T and vector h, required to calculate the corresponding measurement vectors y k , m k , k = 1 , 2 , , K according to (2). Ideally, matrix T would be the 3 × 3 identity matrix while the bias vector h would be the 3 × 1 vector of zeros. A realistic model for T and h, accounting for the sensor’s non-idealities, is derived by using the concept of additive perturbation
T = α I 3 × 3 + E
h = e
where α accounts for gross scaling errors, E is a 3 × 3 perturbation matrix with random, typically small, coefficients and e is 3 × 1 perturbation bias vector with random coefficients. Finally, a sequence of white noise ε N ( 0 , σ 2 ) is added to the measurements and the measurement vectors y k , m k , k = 1 , 2 , , K are derived according to (2)
y = T m + h + ε .
Algorithm 8: Generation of synthetic data
Step 1:Initialize the number of measurements K and the radius of sphere r
Step 2: Calculate Golden Ratio: φ = 1 + 5 2
Step 3: for each  k = 1 , 2 , , K do:
        θ = 2 π k φ
        ϕ = arccos 1 2 ( k 0.5 ) K
        m k = m x , m y , m z = r cos θ sin ϕ , r sin θ sin ϕ , r cos ϕ
Step 4: Pick the scaling parameter, α , the perturbation matrix, E
        and the perturbation vector, e.
Step 5: Calculate T and h according to (66).
Step 6: Generate a sequence of white noise: ε N ( 0 , σ 2 )
Step 7: Calculate the measurement vectors: y k = T m k + h + ε k (2)
The two datasets generated using Algorithm 8 are presented in Figure 1. Note that for visualization purposes, the scaling parameter, α , the perturbation matrix, E, and the perturbation vector, e, used to create each dataset were set to a rather large value.

10.1.1. Experiment Setup and Evaluation Criteria

To evaluate the performance of the algorithms, we used synthetic data, generated by Algorithm 8, and we executed a great number Monte Carlo simulations. Each simulation consisted of 250 runs of each algorithm while in each run, the same dataset was used as input in all algorithms. An uncertainty was introduced in the generation of each dataset by considering a statistical distribution for the elements, E i j , of the perturbation matrix, E, and the elements, e i , of the perturbation vector e (see (66)). Specifically, for the Monte Carlo simulations we assumed
α U [ 0.8 , 1.2 ]
E i j U [ β , β ]
e i U [ γ , γ ]
where β and γ are scalars, the effect of which was tested using multiple Monte Carlo simulations. Note that we considered the scaling factor, α , to be close to the ideal value of α = 1 . That may not be the case when real-world measurements are used, however, it is trivial, and common, to properly scale the measurements before the calibration procedure and remove gross scaling errors. In this way, the algorithms are not burdened, searching for a scaling relationship which can be easily provided by simple data preprocessing.
A challenging point while setting up the experiments was to determine the number of samples of each dataset and the value of the sensor’s noise variance, σ 2 . We considered a dataset of 300 measurements as a solid choice for a simulation environment based on [4,7] while we experimentally confirmed that bigger datasets do not improve the performance of any algorithm. We also examined the performance of the presented algorithms when smaller datasets, consisting of 150 and 50 measurements, are used. As far as the noise variance, σ 2 , is concerned, we considered a nominal value of σ = 0.005 , following [2,4], while we also simulated the cases of more noisy ( σ = 0.05 ) and less noisy ( σ = 0.0005 ) sensors.
The evaluation of the algorithm for each Monte Carlo simulation was completed in terms of accuracy, execution speed, and robustness. We used the execution speed of each algorithm as a metric of computational efficiency and is defined as the inverse of the mean execution time. As a metric of robustness we considered the percentage of datasets for which each algorithm successfully derived a meaningful solution.
The definition of an accuracy metric is a little more involved. Each algorithm was developed to take as inputs the measurement vectors y k , k = 1 , 2 , , K and output the calibration parameters T and h. Comparing the output bias vector h with the true one, h t r u e , which was used in the data generation procedure, was performed by defining the following cost
J h = h t r u e h .
The calibration matrix T on the other hand is derived under a rotational uncertainty and comparing it with the true one, T t r u e , is a more challenging task.
Consider the measurement model of (2). Noting that the true magnetic field vector in (2) is also unknown, and derived by the calibration algorithm, we can write:
y = T t r u e R R T m + h t r u e
where R is an orthogonal matrix in the O ( 3 ) group. Thus, taking into account the rotational invariance of the Euclidean norm which implies that R T m = m , a calibration algorithm may output any matrix T of the form T = T t r u e R . Thus, a proper cost function to compare T and T t r u e is the following
J T = T T t r u e R F
where, the matrix R is defined as the solution of the following minimization problem
R = argmin Ω O ( 3 ) T T t r u e Ω F .
The solution of (72) is given by the orthogonal procrustes problem [33], and it is
R = U V T
where the matrices U and V are derived from the singular value decomposition (SVD) of the matrix T t r u e T T , i.e., T t r u e T T = U Σ V T , where U , V O ( 3 ) and Σ is a diagonal matrix.
Using (69) and (71) we define the following cost function as a metric of accuracy
J = h t r u e h + T T t r u e R F .
Based on the above and given the results of a Monte Carlo simulation consisted of N executions of each algorithm, we define the following metrics of performance:
  • Accuracy is defined as the mean value of the cost J, defined in (74), across all N executions with meaningful output.
  • Mean execution time is defined as the mean value of the execution time of an algorithm.
  • Robustness is defined as the percentage of datasets for which each algorithm successfully derived a meaningful solution.
The robustness criterion can be seen as the frequency in which an algorithm provides a better solution ( T , h ) in the sense of the cost function (74), than the trivial solution ( I 3 × 3 , 0 3 × 1 ) which assumes no bias and non multiplicative errors. Given the cost J o that corresponds to the trivial solution,
J o = h t r u e 0 3 × 1 + I 3 × 3 T t r u e R F
an execution of an algorithm is considered as successful with meaningful output when
J < δ J o
where δ ( 0 , 1 ) is a robustness parameter. If δ is close to 1, it means that only little improvement with respect to J o is sufficient. As δ gets smaller, better solutions are required. Thus, this parameter can be tuned with respect to the test’s objective and the application’s specifications. Given N runs for an algorithm, its robustness is denoted by R B ( % ) and is defined as
R B ( % ) = 1 N i = 1 N U ( J i < δ J o i ) · 100 .
Here J i and J o i are the values of J (74) and J o (75), respectively, corresponding to the ith run of the algorithm and U is a boolean function, which is one if its argument is true and zero otherwise. Let M denote the number of executions meaningful outputs.
Now, the accuracy metric is only applied on the M meaningful outputs according to the robustness test (76), since otherwise the comparison would be unfair for the least stable algorithms. The accuracy of an algorithm over a dataset is denoted by ρ and it is defined as
ρ = 1 M i = 1 N U ( J i < δ J o i ) J i
which is the mean accuracy metric value over the M executions with meaningful outputs.
Similarly, the time-efficiency metric (i.e., mean execution time) is only applied on the M executions with meaningful outputs according to the robustness test (76). Again, this is because otherwise the comparison would be unfair for the least stable algorithms. The mean execution time of an algorithm over a dataset, is denoted by τ and is defined as
τ = 1 M i = 1 N U ( J i < δ J o i ) t i
where t i is the time needed for the i run to be completed. The execution speed of an algorithm is defined as 1 / τ .

10.1.2. Baseline Evaluation

To derive a baseline evaluation of the presented algorithms, we run a Monte Carlo simulation considering typical values for the sensor’s error and noise parameters. In this simulation we neglected the effect of hard-iron and soft-iron distortions which are, in some cases, the dominant terms of the overall error, as well as extreme cases of large manufacturing imperfections. More specifically, 250 different datasets consisting of 300 measurements each, were generated following Algorithm 8 and considering the following distributions of the model disturbances and the measurement noise
α U [ 0.8 , 1.2 ] E i j U [ 0.05 , 0.05 ] e i U [ 0.05 , 0.05 ] σ = 0.005
The distribution ranges in (80) are based on our literature review. The selection β = γ = 0.05 corresponds to the typical case of approximately 5 % distortion for T and bias h. The measurement noise standard deviation is set to a typical value of σ = 0.005 [2,4].
The performance of the seven algorithms is presented in Table 2.

10.1.3. The Effect of the Offset Perturbation Parameter, γ

Under extreme manufacturing imperfections or the effect of hard-iron distortion, the magnitude of the offset vector, h, can be much larger than that in the typical case. In this Section, we examine how larger values of h affect the performance of the presented algorithms. To do so, we run six Monte Carlo simulations, each one comprised of 250 different datasets generated by following Algorithm 8. The offset vector perturbation parameter e i is simulated with gradually increasing magnitude by expanding the selection horizon U [ γ , γ ] . Afterwards, its corresponding impact on each algorithm’s robustness and accuracy is investigated. The distributions of the model disturbances and measurement noise are:
α U [ 0.8 , 1.2 ] E i j U [ 0.05 , 0.05 ] e i U [ γ l , γ l ] σ = 0.005
for various γ
γ = { 0.05 , 0.15 , 0.25 , 0.5 , 0.75 , 1 }
where l = 1 , 2 , , 6 is the index of Monte Carlo simulation. The extreme case of γ = 1 addresses the possibility of bias being clearly comparable and even indistinguishable to the true magnetic vector. Therefore, as γ increases, the algorithms were driven to their limits and their functionality range was identified. All the other parameters were nominal, to ensure a fair comparison. The results of the six Monte Carlo simulations are presented in Figure 2.
The MAG.I.C.AL and the PSO methods are the most robust ones since they function almost always, even for large values of bias, while TWOSTEP and Wu and Shi’s algorithms are a little less stable. In addition, Dorveaux et al. algorithm and EKF seem to be reliable for small to moderate values of bias. All algorithms, except TWOSTEP and EKF are extremely precise when they function properly. No changes in execution speed are noticed, with the exception of MAG.I.C.AL which probably requires more iterations as the bias increases.

10.1.4. The Effect of the Calibration Matrix Perturbation Parameter, β

Similar to the case of the offset vector, h, under extreme manufacturing imperfections or the effect of soft-distortion, matrix T, can also diverge significantly from the typical case of the identity matrix. In this Section, we examine how larger values of perturbation E affect the performance of the presented algorithms. To do so, we run six Monte Carlo simulations, each one based on 250 different datasets generated by following Algorithm 8. The perturbation elements E i j were simulated with gradually increasing magnitude by expanding the distribution range U [ β , β ] . Afterwards, its corresponding impact on each algorithm’s robustness and accuracy is investigated. The distributions of the model disturbances and measurement noise are:
α U [ 0.8 , 1.2 ] E i j U [ β l , β l ] e i U [ 0.05 , 0.05 ] σ = 0.005
for various β
β = { 0.05 , 0.15 , 0.25 , 0.5 , 0.75 , 1 }
where l = 1 , 2 , , 6 is the index of Monte Carlo simulation. As β increases, the algorithms were driven to their limits and their functionality range was identified. All the other parameters were nominal, to ensure a fair comparison. The results of the six Monte Carlo simulations are presented in Figure 3.
The MAG.I.C.AL algorithm and the algorithm of Dorveaux et al. appear to be the most robust and effective, with similar accuracy. The algorithm of Vasconcelos et al., the PSO algorithm and the EKF algorithm succeed only for small to moderate non-orthogonality errors. Vasconcelos et al. achieves accuracy comparable to that of MAG.I.C.AL. The rest of the algorithms tend to fail frequently as these errors increase. What is surprising is that Wu and Shi’s algorithm provides the most accurate solutions for all β values, but with very low robustness. To conclude, most algorithms handle bias distortion better than non-orthogonality errors.

10.1.5. The Effect of Dataset Size, K

In this section, we examine how the dataset size, K, affects the algorithms’ performance. In general, the diversity of the measurement directions is more crucial than the quantity of them, e.g., a dataset of 50 measurements with directions distributed near uniformly on the unit sphere is significantly more suitable for the algorithms than one with thousands of measurements all having approximately the same direction.
According to existing literature [4,5,7], an order of 300 measurements with directions sufficiently covering the unit sphere form an acceptable dataset for the calibration. Here we use datasets with 50, 150, and 300 measurements to test the algorithms’ limits. To do so, we run three Monte Carlo simulations, based on 250 different datasets generated by Algorithm 8. The distributions of the model disturbances and measurement noise are:
α U [ 0.8 , 1.2 ] E i j U [ 0.05 , 0.05 ] e i U [ 0.05 , 0.05 ] σ = 0.005
The dataset size K varied whereas the distributions’ ranges were fixed to nominal to ensure a fair comparison. The results of the three Monte Carlo simulations are presented in Figure 4.
In general, the dataset size, K, does not seem to be important in terms of robustness. Accuracy is surprisingly high even with only 50 measurements, which is probably an outcome of the well distributed measurement directions using the Fibonacci lattice. Furthermore, the algorithms execution time appeared to be linear with K.

10.1.6. The Effect of the Noise Variance, σ 2

In this section, we examine the influence of measurement noise variance σ on algorithms’ robustness and accuracy. The assumption of pure white Gaussian noise in the measurement model was done. We considered a nominal value of σ = 0.005 , following [2,4], while we also simulated the cases of more noisy ( σ = 0.05 ) and less noisy ( σ = 0.0005 ) sensors. With these choices, we represented algorithms’ capabilities under 3 different orders in the magnitude of the error in the measurement. To do so, we run three Monte Carlo simulations, each one based on 250 different datasets generated by following Algorithm 8. The distributions of the model disturbances and measurement noise are:
α U [ 0.8 , 1.2 ] E i j U [ 0.05 , 0.05 ] e i U [ 0.05 , 0.05 ] σ = { 0.0005 , 0.005 , 0.05 }
Finally, all parameters except σ were set to their default ones, to ensure a fair comparison. The results of the three Monte Carlo simulations are presented in Figure 5.
All algorithms appear to be immune to the change of measurement’s output variance σ . What is worth mentioning is that an increase of one order in variance resulted to a decrease of one order in accuracy for most algorithms (i.e., MAG.I.C.AL, Ali et al., Vasconcelos et al., Dorveaux et al., Wu and Shi).

10.2. Algorithm Evaluation Using Real Data

In this section, the aforementioned algorithms are tested using real data. Multiple datasets captured by low-cost magnetic field sensors were used to verify the algorithms’ performance under real-world conditions. In this case, parameters T t r u e and h t r u e are not known in advance. Therefore, the accuracy metric (74) cannot be used. Since, the measurements took place in a specific location, a constant magnitude of magnetic vector, m = 1 was considered. As a result, a proper cost function to evaluate an algorithm’s effectiveness is the following
J r = 1 K i = 1 K m k 2 1 2
where K is the number of measurements and k = 1 , 2 , , K is the measurement index. The estimated magnetic field vector m k for each k is given by
m k = T 1 ( y k h )
where T and h are the outputs of a calibration algorithm. Such a cost function is described by Wu and Shi (52), as well as by Papafotis and Sotiriadis (65).
To evaluate the performance of the presented algorithms, we used two off-the-shelf, low-cost magnetic field sensors, which are typically found in commercial electronic devices, such as smartphones, activity trackers, etc. More specifically, we captured a total of 30 datasets using the LSM9DS1 by STMicroelectronics and the BNO055 by Bosch Sensortec. The operation parameters of the two sensors during the experiment are presented in Table 3.
During the experiment, two sensors were fixed on the same rigid platform which was rotated by hand in several orientations. In Figure 6a, the mean value of the cost function (81) across all the recorded datasets for every algorithm is presented as a metric of accuracy. The robustness of each algorithm, as defined in (77) is presented in Figure 6b. Note that both Figure 6a,b are in agreement with the results obtained in Section 10.1.2 where synthetic data with typical values for sensor’s noise and measurement distortion were considered.

11. Conclusions

To summarize, a complete and extensive study on calibration methods for low-cost magnetometers was carried out by the authors. Seven algorithms were selected for this purpose according to their popularity and their performance. A standard, unified, and complete linear measurement model was used here as the reference model for analyzing all calibration methods. After establishing the full calibration problem, these seven algorithms were discussed and were presented in an easy-to-implement way.
In order to evaluate fairly the presented algorithms’ performance, we proposed a method for optimally generating artificial magnetometer data. This method was used for executing a plethora of Monte Carlo simulations. The evaluation metrics we focused on were the robustness, the accuracy and the efficiency of the algorithms. We designed several experiments to check the behavior of the algorithms under different values in bias, different values in non-orthogonality errors, different number of measurements and finally under various orders of variance in noise. Finally, several datasets of real magnetometer’s data, from two different, low-cost, commercial sensors were used to verify the results obtained using the artificial data.
The following summarizes our findings regarding the studied algorithms and their possible implementation. Except from the objective criteria that we established in Section 10 to evaluate and compare the presented algorithms (accuracy, robustness, computational efficiency), in Table 4 we also evaluate the algorithms in terms of simplicity. Simplicity is used as a (subjective) metric describing our personal experience developing and testing the algorithms. It is related both to the algorithmic complexity of the algorithms (which is not an inherent disadvantage) and the quality of their presentation in the original manuscripts. The algorithms are discussed in chronological order of publication.
TWOSTEP: Extremely time efficient. Works effectively for small distortions. Has low accuracy in general. The method can be generalized to on-orbit calibration.
Crassidis et al.: Easy to implement. Extremely time efficient. Works effectively for small to medium distortions. The method can be generalized to on-orbit calibration. It is the only algorithm that provides online update. It can be considered as a more accurate and effective version of TWOSTEP with similar time complexity.
Dorveaux et al.: Easy to implement. Moderately time efficient. Robust and accurate, but vulnerable to large values of bias.
Vasconcelos et al.: Difficult to implement. Characterized by high time-complexity. Exceptional accuracy and robustness for small distortions.
Ali et al.: Robust and accurate. Very high computational cost. Some prior knowledge of the search space is beneficial. At the beginning of the algorithm, the unknown variables are randomized and, thus, it is not always ensured that the algorithm will reach an optimal point. Thus, a couple of repetitions might be needed. Using modern PSO algorithms which can constrain the search space and handle a few variable inequalities increases the algorithm’s performance significantly.
Wu and Shi: Difficult to implement. Characterized by high time-complexity. Exceptional accuracy even with larger distortion. We noticed a 1% failure of finding an initial estimate due to inadequacy of applying Cholesky decomposition.
MAG.I.C.AL: Easy to implement. Moderately time efficient. Exceptional robustness and accuracy in both synthetic and real data experiments.
To conclude, in this work, we tried to cover a broad range of realistic cases and test the limits of the algorithms, noting that in real life the performance requirements differ from application to another. In some applications computational efficiency may be of major importance while great accuracy may not be needed, while in others, a very accurate calibration is essential even if significantly more computation time is required for this. Thus, there is no “perfect” algorithm appropriate for all applications; different algorithms may be more appropriate for different cases.

Author Contributions

Investigation, K.P., D.N. and P.P.S.; Writing—original draft, K.P., D.N.; Writing—review and editing, K.P., D.N. and P.P.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research is co-financed by Greece and the European Union (European Social Fund- ESF) through the Operational Programme “Human +Resources Development, Education and Lifelong Learning” in the context of the project “Strengthening Human Resources Research Potential via Doctorate Research” (MIS-5000432), implemented by the State Scholarships Foundation (IKY).

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Appendix A.1. Gradient Vector and Hessian Matrix for the Algorithm of Vasconcelos et al.

This section (Section 6) presents the analytical algebraic expressions for the gradient and Hessian of the likelihood function (35), used in descent optimization methods. Let u k = y k h , the gradient of the likelihood function J (33) is denoted by J | x = J | T ^ J | h and described by the submatrices
J | T ^ = k = 0 N 2 c k σ 2 u k T ^ u k
J | h = k = 0 N 2 c k σ 2 T ^ T T ^ u k
where c k = 1 T ^ u k 1 . The Hessian
2 J | x = H T ^ , T ^ H T ^ , h H T ^ , h T H h , h
is given by the following submatrices
H T ^ , T ^ = k = 1 N 2 σ 2 ( u k u k T ) ( T ^ u k u k T T ^ T ) T ^ u k 3 + c k ( u k u k T ) I 3 × 3
H T ^ , h = k = 1 N 2 σ 2 ( u k T ^ u k ) u k T T ^ T T ^ T ^ u k 3 + c k u k T ^ + I 3 × 3 T ^ u k
H h , h = k = 1 N 2 σ 2 T ^ T T ^ u k u k T T ^ T T ^ T ^ u k 3 + c k T ^ T T ^

Appendix A.2. Gradient Vector and Hessian Matrix for the Algorithm of Wu and Shi

This section (Section 8) presents the algebraic expressions for the gradient and Hessian of the likelihood function (51), used in descent optimization methods. For notational simplicity T ^ and m ^ are replaced by T and m. Let u k = y k h , the Jacobian vector and Hessian matrix are, respectively, derived as
J | x = J T T J h T J m k T k = 1 : N J λ k T k = 1 : N T
2 J | x = H T T H T h H T m k 0 9 × 1 H T h T H h h H h m k 0 3 × 1 H T m k T H h m k T H m k m k H m k λ k 0 9 × 1 T 0 3 × 1 T H m k λ k T 0
where
J T = 2 k = 1 N m k u k T m k J h = 2 k = 1 N u k T m k J m k = 2 T T u k T m k + 2 λ k m k J λ k = m k 2 1
and
H T T = 2 k = 1 N m k m k T I H T h = 2 k = 1 N m k I H T m k = 2 m k I T I u k T m k H h h = 2 N I H h m k = 2 T H m k m k = 2 T T T + 2 λ k I H m k λ k = 2 m k
Note that in [4], the calibration matrix, T, is considered to be an upper triangular matrix. Thus, from both the gradient vector and the Hessian matrix, the rows and columns that correspond to the lower triangular elements of T must be removed.

References

  1. Alonso, R.; Shuster, M.D. Complete Linear Attitude-Independent Magnetometer Calibration. J. Astronaut. Sci. 2002, 50, 477–490. [Google Scholar] [CrossRef]
  2. Vasconcelos, J.F.; Elkaim, G.; Silvestre, C.; Oliveira, P.; Cardeira, B. Geometric Approach to Strapdown Magnetometer Calibration in Sensor Frame. IEEE Trans. Aerosp. Electron. Syst. 2011, 47, 1293–1306. [Google Scholar] [CrossRef] [Green Version]
  3. Dorveaux, E.; Vissière, D.; Martin, A.; Petit, N. Iterative calibration method for inertial and magnetic sensors. In Proceedings of the 48h IEEE Conference on Decision and Control (CDC) Held Jointly with 2009 28th Chinese Control Conference, Shanghai, China, 15–18 December 2009; pp. 8296–8303. [Google Scholar] [CrossRef] [Green Version]
  4. Wu, Y.; Shi, W. On Calibration of Three-Axis Magnetometer. IEEE Sens. J. 2015, 15, 6424–6431. [Google Scholar] [CrossRef] [Green Version]
  5. Papafotis, K.; Sotiriadis, P.P. MAG.I.C.AL.—A Unified Methodology for Magnetic and Inertial Sensors Calibration and Alignment. IEEE Sens. J. 2019, 19, 8241–8251. [Google Scholar] [CrossRef]
  6. Crassidis, J.; Lai, K.L.; Herman, R.R. Real-Time Attitude-Independent Three-Axis Magnetometer Calibration. J. Guid. Control Dyn. 2005, 28, 115–120. [Google Scholar] [CrossRef]
  7. Ali, A.; Siddharth, S.; Syed, Z.; El-Sheimy, N. Swarm Optimization-Based Magnetometer Calibration for Personal Handheld Devices. Sensors 2012, 12, 12455–12472. [Google Scholar] [CrossRef] [Green Version]
  8. Papafotis, K.; Sotiriadis, P.P. Accelerometer and Magnetometer Joint Calibration and Axes Alignment. Technologies 2020, 8, 11. [Google Scholar] [CrossRef] [Green Version]
  9. Kok, M.; Hol, J.D.; Schön, T.B.; Gustafsson, F.; Luinge, H. Calibration of a magnetometer in combination with inertial sensors. In Proceedings of the 2012 15th International Conference on Information Fusion, Singapore, 9–12 July 2012; pp. 787–793. [Google Scholar]
  10. Wu, Y.; Zou, D.; Liu, P.; Yu, W. Dynamic Magnetometer Calibration and Alignment to Inertial Sensors by Kalman Filtering. IEEE Trans. Control Syst. Technol. 2018, 26, 716–723. [Google Scholar] [CrossRef] [Green Version]
  11. Kok, M.; Schön, T.B. Maximum likelihood calibration of a magnetometer using inertial sensors. IFAC Proc. Vol. 2014, 47, 92–97. [Google Scholar] [CrossRef] [Green Version]
  12. Li, X.; Li, Z. A new calibration method for tri-axial field sensors in strap-down navigation systems. Meas. Sci. Technol. 2012, 23, 105105. [Google Scholar] [CrossRef]
  13. Cao, G.; Xu, X.; Xu, D. Real-Time Calibration of Magnetometers Using the RLS/ML Algorithm. Sensors 2020, 20, 535. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  14. Hadjigeorgiou, N.; Asimakopoulos, K.; Papafotis, K.; Sotiriadis, P.P. Vector Magnetic Field Sensors: Operating Principles, Calibration and Applications. IEEE Sens. J. 2020, 21, 12531–12544. [Google Scholar] [CrossRef]
  15. IAGA Division V; Working Group 8. Revision of International Geomagnetic Reference Field released. EOS Trans. 1996, 77, 153. [Google Scholar] [CrossRef]
  16. Gambhir, B. Determination of Magnetometer Biases Using Module RESIDG; Technical Report; Computer Sciences Corporation: Falls Church, VA, USA, 1975. [Google Scholar]
  17. LERNER, G. Spacecraft Attitude Determination and Control; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1978. [Google Scholar]
  18. Alonso, R.; Shuster, M.D. TWOSTEP: A fast robust algorithm for attitude-independent magnetometer-bias determination. J. Astronaut. Sci. 2002, 50, 433–451. [Google Scholar] [CrossRef]
  19. Strang, G. Linear Algebra and Its Applications; Brooks Cole/Cengage Learning: Belmont, CA, USA, 2007. [Google Scholar]
  20. Crassidis, J.L.; Markley, F.L.; Lightsey, E.G. Global Positioning System Integer Ambiguity Resolution without Attitude Knowledge. J. Guid. Control Dyn. 1999, 22, 212–218. [Google Scholar] [CrossRef] [Green Version]
  21. Crassidis, J.L. Optimal Estimation of Dynamic Systems; CRC Press: Boca Raton, FL, USA, 2004. [Google Scholar]
  22. Springmann, J.C.; Cutler, J.W. Attitude-Independent Magnetometer Calibration with Time-Varying Bias. J. Guid. Control Dyn. 2012, 35, 1080–1088. [Google Scholar] [CrossRef]
  23. Foster, C.C. Elkaim, Extension of a two-step calibration methodology to include nonorthogonal sensor axes. IEEE Trans. Aerosp. Electron. Syst. 2008, 44, 1070–1078. [Google Scholar] [CrossRef]
  24. Gebre-Egziabher, D.; Elkaim, G.; Powell, J.; Parkinson, B. Calibration of Strapdown Magnetometers in Magnetic Field Domain. J. Aerosp. Eng. 2006, 19, 87–102. [Google Scholar] [CrossRef]
  25. Kennedy, J.; Eberhart, R. Particle swarm optimization. In Proceedings of the ICNN’95—International Conference on Neural Networks, Perth, WA, Australia, 27 November–1 December 1995; Volume 4, pp. 1942–1948. [Google Scholar] [CrossRef]
  26. Kennedy, J.; Obaiahnahatti, B.; Shi, Y. Swarm Intelligence; Morgan Kaufmann Academic Press: San Francisco, CA, USA, 2001; Volume 1, pp. 1931–1938. [Google Scholar]
  27. Magnus Erik Hvass Pedersen. Good Parameters for Particle Swarm Optimization; Hvass Laboratories. 2010. Available online: https://www.semanticscholar.org/paper/Good-Parameters-for-Particle-Swarm-Optimization-Pedersen/a4ad7500b64d70a2ec84bf57cfc2fedfdf770433 (accessed on 31 July 2021).
  28. Mezura-Montes, E.; Coello, C. Constraint-handling in nature-inspired numerical optimization: Past, present and future. Swarm Evol. Comput. 2011, 1, 173–194. [Google Scholar] [CrossRef]
  29. Helwig, S.; Branke, J.; Mostaghim, S. Experimental Analysis of Bound Handling Techniques in Particle Swarm Optimization. IEEE Trans. Evol. Comput. 2013, 17, 259–271. [Google Scholar] [CrossRef]
  30. MATLAB. Optimization Toolbox; The MathWorks Inc.: Natick, MA, USA, 1999. [Google Scholar]
  31. Roberts, M. How to Evenly Distribute Points on a Sphere More Effectively than the Canonical Fibonacci Lattice. Available online: http://extremelearning.com.au/how-to-evenly-distribute-points-on-a-sphere-more-effectively-than-the-canonical-fibonacci-lattice/ (accessed on 22 May 2021).
  32. Gonzalez, A. Measurement of Areas on a Sphere Using Fibonacci and Latitude–Longitude Lattices. Math. Geosci. 2010, 42, 49–64. [Google Scholar] [CrossRef] [Green Version]
  33. Schonemann, P. A generalized solution of the orthogonal procrustes problem. Psychometrika 1966, 31, 1–10. [Google Scholar] [CrossRef]
Figure 1. Two synthetic datasets generated using Algorithm 8 for K = 150 (a) and K = 300 (b), respectively.
Figure 1. Two synthetic datasets generated using Algorithm 8 for K = 150 (a) and K = 300 (b), respectively.
Sensors 21 05288 g001aSensors 21 05288 g001b
Figure 2. Performance characteristics of the presented algorithms for different values of γ .
Figure 2. Performance characteristics of the presented algorithms for different values of γ .
Sensors 21 05288 g002aSensors 21 05288 g002b
Figure 3. Performance characteristics of the presented algorithms for different values of β .
Figure 3. Performance characteristics of the presented algorithms for different values of β .
Sensors 21 05288 g003aSensors 21 05288 g003b
Figure 4. Performance characteristics of the presented algorithms for different values of K.
Figure 4. Performance characteristics of the presented algorithms for different values of K.
Sensors 21 05288 g004
Figure 5. Performance characteristics of the presented algorithms for different values of the noise variance, σ 2 .
Figure 5. Performance characteristics of the presented algorithms for different values of the noise variance, σ 2 .
Sensors 21 05288 g005aSensors 21 05288 g005b
Figure 6. Performance characteristics of the presented algorithms using multiple datasets of real data from two different commercial magnetic field sensors.
Figure 6. Performance characteristics of the presented algorithms using multiple datasets of real data from two different commercial magnetic field sensors.
Sensors 21 05288 g006
Table 1. Notation.
Table 1. Notation.
. Euclidean Norm
. F Frobenius Norm
vec (·)Vectorization of Matrix
diag(·)Diagonal Matrix
chol (·)Cholesky Factorization
I n × n n × n Identity Matrix
0 m × 1 m × 1 Zero Vector
N Normal Distribution
U Uniform Distribution
Gradient Vector
2 Hessian Matrix
Kronecker Product
O ( 3 ) Orthogonal Group of dimension 3
S O ( 3 ) 3D Rotation Group
U ( 3 ) Group of 3 × 3 Upper Triangular Matrices
Table 2. Baseline evaluation of the presented algorithms.
Table 2. Baseline evaluation of the presented algorithms.
AlgorithmAccuracy ( 1 / ρ )Robustness ( RB % ) Execution Speed ( 1 / τ )
TWOSTEP [1] 35.3 × 10 0 91.6 % 455 s 1
Crassidis et al. [6] 3.31 × 10 3 100 % 47.6 s 1
Dorveaux et al. [3] 2.26 × 10 5 100 % 12.8 s 1
Vasconcelos et al. [2] 2.28 × 10 5 99.6 % 0.089 s 1
Ali et al. [7] 2.27 × 10 5 98.8 % 0.10 s 1
Wu and Shi [4] 2.32 × 10 5 87.2 % 0.24 s 1
MAG.I.C.AL [5] 2.28 × 10 5 100 % 29.4 s 1
Table 3. Operation parameters of the two magnetic field sensors.
Table 3. Operation parameters of the two magnetic field sensors.
BNO055LSM9DS1TR
Measurement Range ± 13 Gauss ± 4 Gauss
Sampling Rate30 Hz80 Hz
Measurement Resolution16 bits16 bits
Table 4. Algorithms’ Comparison Summary – More checkmarks correspond to better performance regarding a specific metric.
Table 4. Algorithms’ Comparison Summary – More checkmarks correspond to better performance regarding a specific metric.
AlgorithmSimplicityRobustnessAccuracyEfficiency
TWOSTEP ✓✓ ✓✓ ✓✓✓
Crassidis et al. ✓✓✓ ✓✓ ✓✓✓
Dorveaux et al. ✓✓✓ ✓✓✓ ✓✓✓ ✓✓
Vasconcelos et al. ✓✓
Ali et al. ✓✓ ✓✓✓ ✓✓✓
Wu and Shi ✓✓✓
MAG.I.C.AL ✓✓✓ ✓✓✓ ✓✓✓ ✓✓
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Papafotis, K.; Nikitas, D.; Sotiriadis, P.P. Magnetic Field Sensors’ Calibration: Algorithms’ Overview and Comparison. Sensors 2021, 21, 5288. https://doi.org/10.3390/s21165288

AMA Style

Papafotis K, Nikitas D, Sotiriadis PP. Magnetic Field Sensors’ Calibration: Algorithms’ Overview and Comparison. Sensors. 2021; 21(16):5288. https://doi.org/10.3390/s21165288

Chicago/Turabian Style

Papafotis, Konstantinos, Dimitris Nikitas, and Paul P. Sotiriadis. 2021. "Magnetic Field Sensors’ Calibration: Algorithms’ Overview and Comparison" Sensors 21, no. 16: 5288. https://doi.org/10.3390/s21165288

APA Style

Papafotis, K., Nikitas, D., & Sotiriadis, P. P. (2021). Magnetic Field Sensors’ Calibration: Algorithms’ Overview and Comparison. Sensors, 21(16), 5288. https://doi.org/10.3390/s21165288

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop