IMU calibration

Abstract

A method of calibrating an inertial measurement unit, the method comprising: (a) collecting data from the inertial measurement unit while stationary as a first step; (b) collecting data from the inertial measurement unit while repositioning the inertial measurement unit around three orthogonal axes of the inertial measurement unit as a second step; (c) calibrating a plurality of gyroscopes using the data collected during the first step and the second step; (d) calibrating a plurality of magnetometers using the data collected during the first step and the second step; (e) calibrating a plurality of accelerometers using the data collected during the first step and the second step; (f) where calibrating the plurality of magnetometers includes extracting parameters for distortion detection and using the extracted parameters to determine if magnetic distortion is present within a local field of the inertial measurement unit.

Description

INTRODUCTION TO THE INVENTION

The present disclosure is directed to methods and techniques for calibrating inertial measurement units (IMUs) and, more specifically, to methods and techniques that involve calibrating IMUs to account for local environments within which the IMUs operate.

It should be noted that Patent Cooperation Treaty application PCT/US14/69411, filed Dec. 9, 2014 is hereby incorporated by reference and appended hereto as Appendix A.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a graphical depiction of a ferromagnetic device mounted to an inertial measurement unit and being rotated around an X-axis.

FIG. 2 is a graphical depiction of the ferromagnetic device mounted to the inertial measurement unit of FIG. 1, being rotated around a Y-axis.

FIG. 3 is a graphical depiction of the ferromagnetic device mounted to the inertial measurement unit of FIG. 1, being rotated around a Z-axis.

DETAILED DESCRIPTION

The exemplary embodiments of the present disclosure are described and illustrated below to encompass inertial measurement units and methods of calibrating the same. Of course, it will be apparent to those of ordinary skill in the art that the embodiments discussed below are exemplary in nature and may be reconfigured without departing from the scope and spirit of the present invention. However, for clarity and precision, the exemplary embodiments as discussed below may include optional steps, methods, and features that one of ordinary skill should recognize as not being a requisite to fall within the scope of the present invention.

Many current devices make use of integrated sensor components to capture and track information related to the use, location, and status of the device. One such system, the inertial measurement unit (IMU), is used to track object motion. A common packaging of an IMU comprises a single integrated circuit (IC) consisting of multiple inertial sensors, which can be a combination of accelerometers, gyroscopes, and magnetometers. These three sensor types can also be separated into individual ICs populated on the same circuit board. These two cases (single IC and multiple ICs) can be considered the same for purposes of the instant disclosure.

For all such sensors of the IMU, a calibration step may be performed to ensure the proper referencing and operation in the local environment. For example, accelerometers may be calibrated to determine a reference with respect to the direction of gravity and gyroscopes may be calibrated to some known rotational velocity. Magnetometer calibration is more complex, as the process requires knowledge of the local magnetic field strength and direction at a multitude of magnetometer orientations. The local environment calibration can be carried out by hand or by machine. The minimum considerations for local environment calibration should involve maneuvering of the IMU around the three orthogonal axes of the IMU, and a stationary period.

For applications requiring highly reliable and reproducible calibration of an IMU, a calibration machine may be used. A calibration machine may be constructed so that it does not introduce additional magnetic distortion to the IMUS. The calibration machine should produce movement for calibration for all the sensors of the IMU as well as movement in order for an algorithm to inspect the quality of the calibration.

The process for calibrating measurements from uncalibrated accelerometers, gyroscopes, and magnetometers of an IMU is described hereafter. This exemplary process is applicable to a single IMU or multiple IMUS populated on the same circuit board, where all sensors (uncalibrated accelerometers, gyroscopes, and magnetometers) may be calibrated simultaneously.

In exemplary form, the calibration procedure may include the following routines to allow adequate data to calculate and verify the calibration. This exemplary calibration comprises four stages that include, without limitation: (1) calibrating the gyroscopes; (2) calibrating the magnetometers routine; (3) calibrating the accelerometer routine; and, (4) calibrating the checking/validation routine.

Calibrating the gyroscopes of each IMU should provide for the gyroscopes being stationary for a set amount of time. As part of calibrating the gyroscopes, the instant process is determining the bias of the gyroscope by arithmetic mean as set forth in Equation 1 immediately below:

$\begin{array}{cc}{\mathrm{Gyro}}_{{\mathrm{mean}}_i}\left[\dots \right]={\stackrel{\_}{{\mathrm{DATA}}_{{\mathrm{gyro}}_i}}}_{\mathrm{gyro}\mathrm{calibration}},i=1,2,\dots & \mathrm{Equation}1\end{array}$

Calibrating the magnetometers of each IMU should involve rotating the IMU around all of its orthogonal axes. In this fashion, this magnetometer calibration process first determines the soft and hard iron distortion using multi-dimensional Ellipsoidal Fitting function (EF_fun) set forth in Equation 2 immediately below:

$\begin{array}{cc}{\mathrm{SFe}}_{{\mathrm{orig}}_i},{\mathrm{HFe}}_{{\mathrm{orig}}_i}={\mathrm{EF}}_{\mathrm{fun}\left({\mathrm{DATA}}_{{\mathrm{mag}}_i}\left[\mathrm{Mag}\mathrm{calibration}\mathrm{begin}\to \mathrm{end}\right]\right)},i=1,2,\dots & \mathrm{Equation}2\end{array}$ where: SFe_orig is the transformation of ellipsoid to sphere, and HFe_orig is the sensor and local bias.

The resultant data from Equation 2 are then processed using Equation 3, shown immediately below:

$\begin{array}{cc}{\mathrm{procMAG}}_{{\mathrm{orig}}_i}=\mathrm{processMag}\left({\mathrm{DATA}}_{{\mathrm{mag}}_i}[\mathrm{Acc}\mathrm{calib}r\mathrm{ation}\mathrm{ ration }\mathrm{begin}\to \mathrm{end}\mathrm{of}\mathrm{checking}\mathrm{routine}],{\mathrm{SFe}}_{{\mathrm{orig}}_i},{\mathrm{HFe}}_{{\mathrm{orig}}_i}\right),i=1,2,\dots )& \mathrm{Equation}3\end{array}$

Thereafter, the alignment differences among magnetometers are calibrated. In doing so, Equation 4 may be utilized, with the assumption that all magnetometers are aligned to a reference alignment using this equation. Alternatively, the magnetometers may be aligned to a different reference alignment if desired. In Equation 4, TMAG_ij represents the alignment transformation of magnetometer “j” to reference magnetometer “i”.

$\begin{array}{cc}{\mathrm{TMAG}}_{\mathrm{ij}}={\mathrm{procMAG}}_{{\mathrm{orig}}_j}\backslash {\mathrm{procMAG}}_{{\mathrm{orig}}_i},i=1,j=2,3\dots & \mathrm{Equation}4\end{array}$

Calibrating the accelerometers of each IMU should involve rotating the IMU around two axes perpendicular to gravity. In this fashion, the accelerometer calibration is calculated with the multi-dimensional Ellipsoidal Fitting function (EF_fun) set forth immediately below as Equation 5:

$\begin{array}{cc}{\mathrm{AccScale}}_{{\mathrm{orig}}_i},{\mathrm{AccBias}}_{{\mathrm{orig}}_i}={\mathrm{EF}}_{\mathrm{fun}\left({\mathrm{DATA}}_{{\mathrm{Acc}}_i}\left[\mathrm{Acc}\mathrm{calibration}\mathrm{begin}\to \mathrm{end}\right]\right)},i=1,2,\dots & \mathrm{Equation}5\end{array}$

Calibrating the checking routine should involve rotating the IMU around the axis collinear with gravity. The result is best if the motion of the checking routine is performed multiple times at different angles with respect to gravity.

After calibrating the sensors, one may take into account manufacturing and IC assembly placement errors to extract or account for these errors. In this fashion, the instant disclosure includes a calibration sequence taking into account any mal-alignment among sensors, which are calibrated according to Equation 6-10 (applying the calibration parameters) provided immediately below. In doing so, all sensor data collected during the calibration may be subject to these calibration equations/parameters.

procGyroi = processGyro(DATAgyroi, Gyrobiasi), Equation 6
i = 1, 2, . . .
procAcci = processAcc(DATAAcci, Accbiasi), i = 1, 2, . . . Equation 7
procMAGi = processMag(DATAmagi, SFeorigi, HFeorigi), Equation 8
i = 1, 2, . . .
Continuing the example as aligning all magnetometer to
the first one,
alignProcMAGj = procMAGj * TMAGij, i = 1, j = 2, 3 . . . Equation 9
alignProcMAGi = procMAGi         Equation 10

With multiple IMUs rigidly fixed together, the aligned and calibrated signals should be roughly equal among each sensor type in the case of no noise or isotropic (to the sensor) noise.

Despite calibration with the local environment, magnetometers are especially prone to transient local disturbance. Specifically, introduction of additional ferromagnetic sources or movement of ferromagnetic sources after initial calibration (described previously) can distort the sensor signal and result in inaccurate readings. Distortion states can be classified as hard iron distortion, which is primary caused by the introduction of ferromagnetic materials such as steel tools that can deflect the field lines and produce an offset on the magnetometer sensors, and soft iron distortion, which is caused by an irregular magnetic field (such as magnets) that can alter the shape of the magnetic field. Given that these devices causing distortions will enter and exit the application field (area of use) of the IMU at unpredictable orientations, positions, and times, it is desirable that these distortions be detected and accounted for. By way of example, one area where distortions are commonplace is an operating room during a surgical procedure. More specifically, surgical instruments within an operating room are frequently repositioned within the operating field at different times and locations depending on the stage of the surgical procedure, which necessarily causes magnetic field distortions. It should be noted, however, that the exemplary disclosure is not limited to calibration of IMUs for surgical navigation or any particular field of use. Rather, the instant disclosure is applicable to calibrating any IMU in any field of use where local distortions may be present.

There are two essential parts to detect transient magnetic distortion. The first part is that the IMU must comprise two or more magnetometers (that may be populated on a common circuit board). The magnetometers are calibrated and aligned based the foregoing calibration teachings. The second part is the detection algorithm, which consists of two processes: (1) extracting parameters for distortion detection, which are calculated from the calibration data; and, (2) use of extracted parameters on captured IMU data to check if distortion is present.

What follows is a detailed explanation for calculating appropriate distortion threshold parameters. In general, one may use an algorithm for this calculation that may make use of the following assumptions: (1) IMU motion (or lack of motion) is detectable by all IMUs (accelerometer, gyroscope and magnetometer); (2) the angle and length of the calibrated magnetometers should have approximately identical values; (3) the calibrated magnetometer vectors should be approximately unity in magnitude; (4) the angle between the calibrated magnetometer vector and calibrated accelerometer vector should not change; (5) the radius of the calibrated magnetometers should be approximately equal to one; and, (6) the quaternions calculated using different magnetometers should be approximately equal.

Sufficient deviation from assumption 1 is deemed a motion anomaly and deviations from assumptions 2-6 are deemed a magnetic anomaly. Parameters for motion detection are calculated for each IMU sensor (accelerometer, gyroscope, and magnetometer) using data collected during local environment calibration. Using stationary and dynamic data, appropriate thresholds of motion are calculated for every IMU sensor, so that the calibrated IMU sensor output above the threshold indicates that the sensor is in motion and below which the sensor can be considered stationary (see Equation 2). A motion anomaly is considered to occur as in Table (2), which is described in more detail hereafter. In short, if a magnetometer has detected motion, but other sensors do not, a motion anomaly is said to exist.

For the magnetic anomaly, a deviation from assumption 2 is used to create a cost function, which is calculated as the magnitude of the vector difference weighted with the angular derivation between the vectors (see Algorithm 1 and Algorithm 4). During local environment calibration, expected values for the deviation from unity and from identical signals are calculated to form a cost function. The threshold is calculated from the value of this cost function during local environment calibration (when no additional distortion is presented). Then, as part of an anomaly detection process, any values of this cost function above the threshold are considered a magnetic anomaly.

For assumption 3, the magnetic strength is calculated as the mean of combined magnitudes (Algorithm (5)), which without noise is expected to be approximately 1. Then, as part of the anomaly detection process, any value outside the calculated magnetic strength (including a preset tolerance) is considered an anomaly.

For assumption 4, the derivation is calculated as the angular derivation between the magnetometer and the accelerometer vectors (Algorithm (6)). Then, as part of the anomaly detection process, any value outside the calculated angle (including a preset tolerance) is considered an anomaly.

For assumption 5, the deviation is calculated by using all of the magnetometer readouts, along with a preset of artificial points, to estimate a radius of the magnetic field using an ellipsoid estimation function as shown in Algorithm (7). The artificial points are used to support the estimation calculation. However, if the IMU contains more than 11 magnetometers, no artificial point is needed to resolve the radius. Then, as part of the anomaly detection process, any value outside the calculated radii (including a preset tolerance) is considered an anomaly.

For assumption 6, the deviation is determined by first calculating the orientation of the IMU using sensor fusion algorithm (e.g. Extended Kalman Filter) using different magnetometers as input, and an angular difference among the output orientation is calculated as shown in Algorithm (8). Then, as part of the anomaly detection process, any values greater than the angular different with a preset tolerance is considered an anomaly.

If motion or magnetic anomaly occurs, some unaccounted for magnetic distortion is present. An alert system for distortion provides feedback indicating the reliability of the orientation reported by the sensors may be inaccurate. What follows is a more detailed discussion of the algorithms utilized to detect motion and magnetic anomalies.

The anomaly detection is based on three categories of detection. The first is sensor fault detection, the second is magnetometer radii verification, and the third is quaternion output discrepancies. Sensor fault detection consists of the following functions: (i) motion and magnetic anomalies; (ii) magnetic strength; and (iii) angle between magnetometer and accelerometer vectors. The following method is used to consistently extract threshold values (TVs) for motion anomalies and cost functions from the calibration. The method utilizes Algorithm (1) as follows:

Algorithm (1)
Iteration: for [beginning of gyro calibration → end of checking routines]
Gyromeani[. . .] = |procGyrocurrent-kernelcurrent|, i = 1, 2, . . .
Gyrostdi[. . .] = std(procGyro[current − kernel → current]), i = 1, 2, . . .
Accstdi[. . .] = std(procAcci[current − kernel → current]), i = 1, 2, . . .
Magstdi[. . .] = std(alignProcMAGi[current − kernel → current]), i = 1, 2, . . .
lenDiffi,j[⁢...⁢]=∑axis⁢(∑current-kernelcurrent⁢(alignProcMAGi-alignProcMAGj)2)kernel,
i ≠ j, i = 1, 2 . . . , j = 1, 2 . . .
angDiff                                                      i                            ,                            j                                                                          [                                                                                                  ⁢                                                  .                                                                                                          .                                                                                                          .                                                ⁢                                                                                                  ]                                            =                                                                                                                              1                            -                                                          (                                                                                                                                    alignProcMAG                                    _                                                                                                        ι                                                                          current                                      -                                      kernel                                                                                                        current                                                                ·                                                                                                      alignProcMAG                                    _                                                                                                        j                                                                          current                                      -                                      kernel                                                                                                        current                                                                                            )                                                                                                                                                                                      norm                            ⁡                                                          (                                                                                                alignProcMAG                                  _                                                                                                  ι                                                                      current                                    -                                    kernel                                                                                                  current                                                            )                                                                                ·                                                      norm                            ⁡                                                          (                                                                                                alignProcMAG                                  _                                                                                                  j                                                                      current                                    -                                    kernel                                                                                                  current                                                            )                                                                                                                                            ,
i ≠ j, i = 1, 2 . . . , j = 1, 2 . . .
metricLenAngi,j[. . .] = lenDiffi,j[. . .] * angDiffi,j[. . .], i ≠ j, i = 1, 2 . . . , j = 1, 2 . . .
currentMagVal                        ⁢                                                                                                  [                                                                                                  ⁢                                                  .                                                                                                          .                                                                                                          .                                                ⁢                                                                                                  ]                                            =                                              1                                                                              max                            ⁡                                                          (                              i                              )                                                                                ⁢                                                      (                                                          ∑                                                                                                                                    ∑                                                                          current                                      -                                      kernel                                                                        current                                                                    ⁢                                                                      (                                                                                                                  alignProcMAG                                        ι                                                                            *                                                                              alignProcMAG                                        j                                                                                                              )                                                                                                                                                        )                                                                                                                ,
i ≠ j, i = 1, 2 . . . , j = 1, 2 . . .
End

Once the iteration from Algorithm (1) is completed, the mu and sigma of the cost function can be determined using Equations 11 and 12, immediately recited below.

μ = currentMagValmag and acc calibration Equation 11
δ = SSF1 * std(currentMagVal)mag and acc calibration Equation 12
Where: (a) SSF1 is a sensitivity scale factor; and,
(b) δ serves as the detection tolerance based on standard
derivation (i.e., how many standard deviations away
from the mean is considered acceptable).

Algorithm (2) is utilized to calculate motion and no motion (stationary) states for the sensors.

Algorithm (2)
Stationarygyrostd  = Gyrostdgyro calibration
Motiongyrostd  = Gyrostdacc calibration
Stationarygyromean  = Gyromeangyro calibration
Motiongyromean  = Gyromeanacc calibration
Stationary                                                      Acc                            i                                                                          =                                                                                                            Acc                                                              ι                                std                                                                                      _                                                                                                                                                                                                    gyro                              ⁢                                                                                                                          ⁢                              calibration                                                                                                                          ,                                              i                        =                        1                                            ,                      2                      ,                                              .                                                                                                  .                                                                                                  .                                                              ⁢
Motion                                                      Acc                            i                                                                          =                                                                                                            Acc                                                              ι                                std                                                                                      _                                                                                                                                                                                                    acc                              ⁢                                                                                                                          ⁢                              calibration                                                                                                                          ,                                              i                        =                        1                                            ,                      2                      ,                                              .                                                                                                  .                                                                                                  .                                                              ⁢
Stationary                                                      MAG                            i                                                                          =                                                                                                            MAG                                                              ι                                std                                                                                      _                                                                                                                                                                                                    gyro                              ⁢                                                                                                                          ⁢                              calibration                                                                                                                          ,                                              i                        =                        1                                            ,                      2                      ,                                              .                                                                                                  .                                                                                                  .                                                              ⁢
Motion                                                      MAG                            i                                                                          =                                                                                                            MAG                                                              ι                                std                                                                                      _                                                                                                                                                                                                    Checking                              ⁢                                                                                                                          ⁢                              routine                                                                                                                          ,                                              i                        =                        1                                            ,                      2                      ,                                              .                                                                                                  .                                                                                                  .                                                              ⁢

A check is performed here to ensure all motion states should have values greater than the corresponding stationary states.

Algorithm (3) is utilized to calculate TVs for motion anomalies.

Algorithm (3)
TVgyrostd = SSF2 * (Motiongyrostd − Stationarygyrostd) + Stationarygyrostd
TVgyromean = SSF3 * (Motiongyromean − Stationarygyrostd) + Stationarygyromean
TVAcci = SSF4 * (MotionAcci − StationaryAcci) + StationaryAcci,
i = 1, 2, . . .
TVMAGi = SSF5 * (MotionMAGi − StationaryMAGi) + StationaryMAGi,
i = 1, 2, . . .
SSF2 to SSF5 are the sensitivity scale factor. It is typically set between 0 and 1, and should be adjusted based on sensors and desired sensitivity of the detection method.

Determination of a magnetic anomaly may be based upon a cost function TV, which is accomplished by calculating the cost function over the period of the calibration period as set forth in Algorithm (4).

Algorithm (4)
Iteration: for [beginning of gyro calibration → beginning of Check 1]
pdfVal⁢[⁢...⁢]=12⁢⁢π*δ⁢e-0.5*((1-μ)δ)2-12⁢⁢π*δ⁢e-0.5*((currentMagVal⁢[⁢...⁢]-μ)δ)2
costFunction⁢[⁢...⁢]=SSF⁢⁢6*pdfVal⁢[⁢...⁢]*∑metricLenAngi,j[⁢...⁢]max⁡(i),i≠j,⁢i=1,2⁢...⁢,j=1,2⁢...⁢
end
TVcostFunction  = max(costFunction) + (costFunction) + std(costFunction)
SSF6 is a scale factor that is used to prevent losing the precision of the cost function from the
small number.

In order to consistently extract TVs for discerning the strength of a magnetic anomaly from the calibration data, one can use Algorithm (5).

Algorithm (5)
Iteration: for [beginning of mag calibration → end of mag calibration]
MagStrii⁢[⁢...⁢]=∑axis⁢(alignProcMAGi)⁢⁢i=1,2,...⁢
end
min_MagStri = min(MagStri) i = 1, 2, . . .
max_MagStri = max(MagStri) i = 1, 2, . . .

In order to consistently extract TVs for discerning the angle between magnetometer and accelerometer vectors from the calibration data, one can use Algorithm (6).

Algorithm (6)
Iteration: for [beginning of gyro calibration → end of checking routine]
Use Extended Kalman Filter to estimation the orientation during the
calibration
qi[. . . ] = EKF(procGyro, procAcc, alignProcMAGi), i = 1, 2, . . .
ei[ . . . ] = [q2i  q3i  q4i], i = 1, 2, . . .
e^i[⁢...⁢]=[0-q⁢⁢4iq⁢⁢3iq⁢⁢4i0-q⁢⁢2i-q⁢⁢3iq⁢⁢2i0],i=1,2,...⁢
Rpi[. . . ] = (q1i2 − ei′ei) * I3×3 + 2(eiei′) − 2q1iêi, i = 1, 2, . . .
MFi[. . . ] = Rpi * alignProcMAGi, i = 1, 2, . . .
GFi[. . . ] = Rpi * proAcc, i = 1, 2, . . .
∠GMi[⁢...⁢]=a⁢cos⁡[MFi·GFiMFi*GFi],i=1,2,...
end
min∠TV_GMi = ∠GMι − tol1
max∠TV_GMi = ∠GMι + tol2
Where: to/1 and to/2 are tolerances for the detection threshold.

A magnetometer radii check calculation is performed on the calibrated magnetometer data before alignment. The radii of a normal point set should be equal to one. Algorithm (7) may be used to determine the radii threshold. As part of using Algorithm (7), a set of artificial points on a unit sphere are created or generated, where the set of points should be close to evenly distributed and have a minimum of 18 points. Using an ellipsoid fitting method on the artificial points, one can utilize the processed magnetometer data to determine the radii as set forth below

Algorithm (7)
Radii = ellipsoidFit([artificial points, procMAGi], i = 1, 2 . . .
minRadii = Radii − tol3
maxRadii = Radii + tol4
Where: tol3 and tol4 are tolerances for the detection threshold; and,
minRadii and maxRadii should be within 0.93 and 1.07 for
applications that require high reliability.

Quaternion output discrepancies is a detection method that is based on the orientation output calculated from using different magnetometers. The following algorithm, Algorithm (8), may be used to determine the quaternion output discrepancies threshold, using the quaternions output from Algorithm (6), and calculate the angle between the quaternions.

Algorithm (8)
qDiffi, j[ . . . ] = quatAngleDiff(qi[ . . . ], qj[ . . . ]), i ≠ j, i = 1, 2 . . . ,
j = 1, 2 . . .
TV_qDiffi, j = qDiffi, j[ . . . ] + SSF7 * std(qDiffi, j[ . . . ]), i ≠ j,
i = 1, 2 . . . , j = 1, 2 . . .

After extracting the TVs for anomaly detection, the following calculation may be used to determine if the IMU is being affected by a local anomaly. As a prefatory step, all incoming sensor data may be processed using Equations 1-10 discussed previously. In addition, Algorithms (1) and (4) may be utilized to determine relative deviation among magnetometers within a preset kernel size. For the exemplary sensor fault detection system, the motion is detected based on the following rules in Table (1).

TABLE(1)
Motion Detection based on sensors
Activity	Condition	Result
Gyro has detected motion	$\mathrm{If}\mathrm{std}\left({\mathrm{procGyro}}_i\right)>{\mathrm{TV}}_{{\mathrm{gyro}}_{{\mathrm{std}}_i}},i=1,2,...$	inertial_motion_detection = 1
	$\mathrm{OR}{\stackrel{\_}{\mathrm{procGyro}}}_i>{\mathrm{TV}}_{{\mathrm{gyro}}_{{\mathrm{mean}}_i}},i=1,2,...$
Accelerometer has	If std(procAcci) > TVAcci , = 1, 2, . . .	inertial_motion_detection = 1
detected motion
Magnetometer has	If std(procMAGi) > TVMAGi , i = 1, 2, . . .	Magnetic_motion_detection = 1
detected motion

Within the preset kernel size, the following rules in Table (2) are used to detect a motion anomaly.

TABLE (2)
Motion anomaly detection
Motion Anomaly Flag Condition
False               If Σinertial_motion_detection > 1
                    AND
                    Σmagnetic_motion_detection > 1
False               If Σinertial_motion_detection < 1
                    AND
                    Σmagnetic_motion_detection < 1
True                If Σinertial_motion_detection < 1
                    AND
                    Σmagnetic_motion_detection > 1
False               If Σinertial_motion_detection > 1
                    AND
                    Σmagnetic_motion_detection < 1

The following rules in Table (3) are used to detect a magnetic anomaly.

TABLE (3)
Magnetic Anomaly detection
Magnetic Anomaly Flag Condition
True                  costFunction > TVcostFunction
False                 costFunction ≤ TVcostFunction

In accordance with the instant disclosure, one can calculate the magnetic strength using Algorithm (5). Thereafter, using the rules of Table (4), one can discern whether a magnetic anomaly is present.

TABLE (4)
Magnetic strength anomaly detection
Magnetic
Strength
Anomaly
Flag  Condition
True  MagStri < min_MagStri OR MagStri > max_MagStri
False MagStri ≥ min_MagStri OR MagStri ≤ max_MagStri

In accordance with the instant disclosure, one can determine the angle between magnetometer and accelerometer vectors using Algorithm (6). This detection is only activated when the IMU/sensors is/are stationary, which is detected from the motion detection algorithm. Using the following table, Table (5), one can discern whether a magnetometer and accelerometer vector anomaly is present.

TABLE (5)
magnetometer and accelerometer vectors anomaly detection
MAG and
ACC vectors
Anomaly Flag Condition
True         LGMi < min   TV_GMi OR   GMi >
             max   TV_GMi
False        LGMi ≥ min   TV_GMi OR   GMi ≤
             max   TV_GMi

In accordance with the instant disclosure, one can determine the magnetic radii by processing the raw data from the magnetometers using Algorithm (7). Using the following table, Table (6), one can discern whether a magnetic radii anomaly is present.

TABLE (6)
Magnetic Radii anomaly detection
Magnetic Radii
Anomaly Flag Condition
True         radii < minRadii OR radii > maxRadii
False        radii ≥ minRadii OR radii ≤ maxRadii

In accordance with the instant disclosure, one can calculate the orientation of the IMUs with Algorithm (6) and calculate the angular difference with Algorithm (8). Using the following table, Table (7), one can discern whether a quaternion output discrepancy is present.

TABLE (7)
Quaternion discrepancies detection
Quaternion
Discrepancies Flag Condition
True               qDiffi,j > TV_qDiffi,j
False              qDiffi,j ≤ TV_qDiffi,j

Anomaly detection can be weighted using the output from each of these sub-detection systems, and determine the relevant level of disturbance detected by the system. For systems requiring high levels of reliability, it should be considered an anomaly if any of these subsystem flags is detected.

The following discussion describes a process for calibrating an IMU when it is being used with or in proximity to a ferromagnetic object. One example, in the medical device field of use, is during surgical navigation, where the IMU is needed to track the motion of a ferromagnetic object (e.g., tools made from stainless steel, CoCr, etc). Since local environment calibration is performed presumably without these ferromagnetic objects in proximity, there is no way to compensate for their distortions during the local environment calibration stage. A secondary calibration step should be undertaken to correct the distortion when ferromagnetic objects will be used in proximity to the IMUs. The proposed method accounts for the fact that an IMU may be rigidly fixed to the ferromagnetic object when in use.

While it is possible to manually maneuver the IMU attached to the ferromagnetic object to collect a lump sum of data for re-calibration, this method is inefficient and requires a substantial amount of time for the user to repeat the maneuver for each ferromagnetic object needing calibration. Secondly, it is possible that the calibration of the combined effect of the local environment and the ferromagnetic object is not aligned with the calibration of the local environment alone. Though either of the calibrations may be functional for yaw tracking, the latent transformation between calibrations on the reference field can produce errors. By way of example, one may think of it as comparing the reading of two compasses, one being horizontal to the ground, the other being tilted 5 degrees. The horizontal compass is calibrated to the local environment, while the other compass is calibrated with both ferromagnetic objects and local environment. While both calibrations are both functional, the reading from each compass will be different. Accordingly, an exemplary algorithm of the instant disclosure addresses both of these issues.

The instant approach makes use of the distortion characteristic of low field strength ferromagnetic objects, and models the correction based on a limited amount of data input. This enables a much simpler calibration maneuver, while still being able to achieve reasonable calibration.

Referring to FIGS. 1-3, the object 100 attached to the IMU 102 (or the IMU 102 by itself) is moved through a series of orientations—which can be confirmed from the IMU data. From the collected data, an initial estimation of the distortion field can be calculated through standard ellipsoid processing and analysis. A reference dataset (currently using 25000 vectors) is drawn from a von-Mises Fisher density, thereby creating a set of uniformly distributed vectors around the IMU sphere. The magnetic distortion calculated during the local environment calibration, and the initial estimation of magnetic distortion with the ferromagnetic object are applied in reverse to the reference dataset. This creates two datasets: (1) a reference dataset being distorted by the local environment; and, (2) a reference dataset being distorted by both the local environment and the attached ferromagnetic object 100 (or ferromagnetic object in proximity to the IMU 102). Since all the points between these two datasets originate from the same reference dataset, a point-to-point correspondence between local environment distortion, and local environment distortion with the ferromagnetic object is established. The correction for the ferromagnetic object is then extracted from the latter using a recursive optimization algorithm with the prior dataset.

This exemplary calibration process can be used for any number of objects, so that multiple objects which have latent magnetic fields can be used with the magnetometer without compromising sensor accuracy, as well as maintain a consistent reference field.

The following section outlines the algorithm for calibrating a ferromagnetic object with a pre-calibrated IMU. This part assumes that the IMU has been calibrated according to the foregoing discussion. The IMU, as part of this calibration sequence, should be rigidly attached to the tool or otherwise have the tool in a constant position with respect to the IMU. The operator rotates the object for at least one revolution along each axis in a standard Cartesian system as depicted in FIGS. 1-3. Once the tool calibration dataset is collected, the following algorithm can calculate and compensate the distortion.

As an initial matter, an initial guess is made as to the center offset of the data set using Equation 13 below.

$\begin{array}{cc}{\mathrm{SFe}}_{{\mathrm{inst}}_i},{\mathrm{HFe}}_{{\mathrm{inst}}_i}={\mathrm{EF}}_{\mathrm{fun}\left({\mathrm{DATA}}_{{\mathrm{mag}}_i}\right)},i=1,2,\dots & \mathrm{Equation}13\end{array}$

Thereafter, a large number (e.g., 25,000) of random 3D unit vectors are generated and application of the inverse of the original and instrument SFe on to them using Algorithm (9) below.

Algorithm (9)
Iteration: for 1 . . . 25000
samples = rand3Dvectors
samplesdistortedi = samples * inv(SFeinsti), i = 1, 2, . . .
samplesnormali = samples * inv(SFeorigi), i = 1, 2, . . .
End

Using Equation 14, shown below, calculate the point correspondence transformation.T_iD→N=samples_normali\samples_distortedi, i=1,2,  Equation 14

Using Equation 15, shown below, calculate the point-to-point distortion distance.Diff_NDi=samples_normali−samples_distortedi, i=1,2,  Equation 15

Using Equation 16, shown below, for each axis, locate the maximum negative of Diff_ND.loc_i=minmax(Diff_NDi), i=1,2,  Equation 16

Using Equation 17, shown below, calculate the index location of loc.

The index location of loc isind_i=index (minmax(Diff_NDi)), i=1,2,  Equation 17

Using Equation 18, shown below, calculate the center along the maximum distortion distance axes.c_i=min(loc_i)+max(loc_i), i=1,2,  Equation 18

Using Equation 19, shown below, calculate the maximum correction boundaries.Dbound_i=samples_distortedi(ind_i), i=1,2,  Equation 19

Using Equation 20, shown below, after the distortion compensation parameters are calculated, center the instrument calibration dataset.centerMAG_i=DATA_magi−HFe_insti, i=1,2,  Equation 20

The foregoing calibration correction can be applied in at least two ways. In a first fashion, the calibration correction is applied using a point correspondence transformation, such as Equation 21 shown below.compensatedMAG_i=centerMAG_i*inv(T_iD→N)i=1,2,  Equation 21

Alternatively, in a second fashion, the calibration correction is applied using geometric scaling and compensation, such as by using Equations 22 and 23.

$\begin{array}{cc}{\mathrm{SF}}_i\left(X,Y,Z\right)=\frac{\begin{array}{c}{\mathrm{centerMAG}}_i\left(X,Y,Z\right)-c\left(X,Y,Z\right)·\\ {\mathrm{Dbound}}_i\left(X,Y,Z\right)-c\left(X,Y,Z\right)\end{array}}{\begin{array}{c}\mathrm{norm}\left({\mathrm{centerMAG}}_i\left(X,Y,Z\right)-c\left(X,Y,Z\right)\right)·\\ \mathrm{norm}\left({\mathrm{Dbound}}_i\left(X,Y,Z\right)-c\left(X,Y,Z\right)\right)\end{array}},i=1,2,\dots & \mathrm{Equation}22\\ {\mathrm{ccompensatedMAG}}_i={\mathrm{centerMAG}}_i+{\mathrm{loc}}_i*{\mathrm{SF}}_i,i=1,2,\dots & \mathrm{Equation}23\end{array}$

Using equation 24, the adjusted calibration parameters may be calculated.SFe_adjInsti, HFe_adjIinsti=magCal(compensatedMAG_i), i=1,2,  Equation 24

Post conclusion of the foregoing calibration computations, the calibration computations may be applied to the magnetometers of each IMU that is attached to the ferromagnetic object or will be in a fixed positional relationship thereto. In order to apply the calibration computations to the magnetometers of each IMU, the magnetometer data is centered using Equation 20. By way of example, on may apply the correction to account for the distortion cause by the tool via a point correspondence transformation using Equation 21, or apply the geometric scaling and compensation using Equations 22 and 23. Thereafter, the magnetometers of the IMU can be processed by applying the adjusted soft and hard iron compensation parameters via Equation 25, shown below.procMAG_i=processMag(compensatedMAG_i, SFe_adjInsti, HFe_adjIinsti), i=1,2,  Equation 25

Once the magnetometers have been processed, the magnetometers are aligned using Equations 9 and 10

Following from the above description, it should be apparent to those of ordinary skill in the art that, while the methods and apparatuses herein described constitute exemplary embodiments of the present invention, the invention described herein is not limited to any precise embodiment and that changes may be made to such embodiments without departing from the scope of the invention as defined by the claims. Additionally, it is to be understood that the invention is defined by the claims and it is not intended that any limitations or elements describing the exemplary embodiments set forth herein are to be incorporated into the interpretation of any claim element unless such limitation or element is explicitly stated. Likewise, it is to be understood that it is not necessary to meet any or all of the identified advantages or objects of the invention disclosed herein in order to fall within the scope of any claims, since the invention is defined by the claims and since inherent and/or unforeseen advantages of the present invention may exist even though they may not have been explicitly discussed herein.

Claims

1. A method of calibrating an inertial measurement unit, the inertial measurement unit comprising at least two magnetometers, and at least two accelerometers, the method comprising: presuming at least one of the following: (1) a vector angle and a vector length of each of the at least two magnetometers have identical values; (2) each of a plurality of magnetometer vectors is one in magnitude; (3) an angle between a calibrated magnetometer vector and a calibrated accelerometer vector should not change; (4) a radius of the at least two magnetometers, post calibration, is equal to one; and (5) quaternions calculated using different ones of the at least two magnetometers are equal;collecting data from the inertial measurement unit, while stationary as a first step;collecting data from the inertial measurement unit while repositioning the inertial measurement unit as a second step;calibrating the at least two magnetometers using the data collected during the first step and the second step; andcalibrating of the at least two accelerometers using the data collected during the first step and the second step.

2. The method of claim 1, wherein at least two of the following are presumed: (1) a vector angle and a vector length of each of the at least two magnetometers have identical values; (2) each of a plurality of magnetometer vectors is one in magnitude; (3) an angle between a calibrated magnetometer vector and a calibrated accelerometer vector should not change; (4) a radius of the at least two magnetometers, post calibration, is equal to one; and (5) quaternions calculated using different ones of the at least two magnetometers are equal.

3. The method of claim 1, wherein at least three of the following are presumed: (1) a vector angle and a vector length of each of the at least two magnetometers have identical values; (2) each of a plurality of magnetometer vectors is one in magnitude; (3) an angle between a calibrated magnetometer vector and a calibrated accelerometer vector should not change; (4) a radius of the at least two magnetometers, post calibration, is equal to one; and (5) quaternions calculated using different ones of the at least two magnetometers are equal.

4. The method of claim 1, wherein at least four of the following are presumed: (1) a vector angle and a vector length of each of the at least two magnetometers have identical values; (2) each of a plurality of magnetometer vectors is one in magnitude; (3) an angle between a calibrated magnetometer vector and a calibrated accelerometer vector should not change; (4) a radius of the at least two magnetometers, post calibration, is equal to one; and (5) quaternions calculated using different ones of the at least two magnetometers are equal.

5. The method of claim 1, wherein all five of the following are presumed: (1) a vector angle and a vector length of each of the at least two magnetometers have identical values; (2) each of a plurality of magnetometer vectors is one in magnitude; (3) an angle between a calibrated magnetometer vector and a calibrated accelerometer vector should not change; (4) a radius of the at least two magnetometers, post calibration, is equal to one; and (5) quaternions calculated using different ones of the at least two magnetometers are equal.

6. The method of any one of claims 1-5, wherein presuming the vector angle and the vector length of the at least two magnetometers have identical values, the method further comprising: calculating a magnitude of vector angle differences weighted using angular derivation between the vector angles to form a cost function;calculating a threshold using the cost function; andestablishing a presence of a magnetic anomaly based upon the threshold.

7. The method of any one of claims 1-5, wherein presuming each of the plurality of magnetometer vectors is one in magnitude, the method further comprising: calculating magnetic strength as a mean of combined magnetic strengths from the at least two magnetometers;establishing a tolerance for the calculated magnetic strength to create a magnetic strength range; andestablishing a presence of a magnetic anomaly based upon the magnetic strength range.

8. The method of any one of claims 1-5, wherein presuming the angle between the calibrated magnetometer vector and the calibrated accelerometer vector should not change, the method further comprising: calculating angular derivation between a magnetometer vector and a accelerometer vector;establishing a tolerance for the calculated angular derivation to create an angular derivation range; andestablishing a presence of a magnetic anomaly based upon the angular derivation range.

9. The method of any one of claims 1-5, wherein presuming the radius of the at least two magnetometers, post calibration, is equal to one, the method further comprising: estimating a radius of a magnetic field using an ellipsoid estimation function and data from the at least two magnetometers;establishing a tolerance for the estimated radius of the magnetic field to create a magnetic field radii range; andestablishing a presence of a magnetic anomaly based upon the magnetic field radii range.

10. The method of any one of claims 1-5, wherein presuming quaternions calculated using different ones of the at least two magnetometers are equal, the method further comprising: calculating an orientation of the inertial measurement unit in quaternion for each of the at least two magnetometers using outputs from the at least two magnetometers;determining a quaternion output discrepancies threshold using the calculated orientations of the inertial measurement unit in quaternion;calculating an angle between the calculated orientations of the inertial measurement unit in quaternion;establishing a tolerance for the calculated angle to create a quaternion angle range; andestablishing a presence of a magnetic anomaly based upon the quaternion angle range.

11. A method of calibrating at least two magnetometers of an inertial measurement unit, while the inertial measurement unit is near an operating room of a medical facility comprising a local environment, the inertial measurement unit including at least two accelerometers, the method comprising: collecting data from the at least two magnetometers of the inertial measurement unit while the inertial measurement unit is rigidly attached to a ferromagnetic object and the ferromagnetic object is repositioned in the local environment;calculating magnetic distortion from the data collected from the at least two magnetometers during repositioning the ferromagnetic object in the local environment;applying the magnetic distortion to a reference dataset, the reference dataset comprising uniformly distributed vectors around a sphere, to create a first distorted dataset;applying an initial magnetic distortion to the reference dataset, the initial magnetic distortion being from a prior calibration of the at least two magnetometers when not rigidly attached to the ferromagnetic object, to create a second distorted dataset; andcalibrating the at least two magnetometers to account for ferromagnetic distortion using a recursive optimization algorithm applied to both the first and second distorted datasets to extract the ferromagnetic distortion.

12. The method of claim 11, wherein the reference dataset comprises thousands of random three dimensional unit vectors drawn from a vin-Mises Fisher distribution around a sphere; and wherein application of the initial magnetic distortion and the magnetic distortion are applied in reverse to the reference dataset.

13. The method of claim 11, wherein the prior calibration of the at least two magnetometers includes: presuming at least one of the following: (1) a vector angle and a vector length of each of the at least two magnetometers have identical values; (2) each of a plurality of magnetometer vectors is one in magnitude; (3) an angle between a calibrated magnetometer vector and a calibrated accelerometer vector should not change; (4) a radius of the at least two magnetometers, post calibration, is equal to one; and (5) quaternions calculated using different ones of the at least two magnetometers are equal;collecting data from the inertial measurement unit, while stationary as a first step;collecting data from the inertial measurement unit while repositioning the inertial measurement unit as a second step;calibrating the at least two magnetometers using the data collected during the first step and the second step; andcalibrating of the at least two accelerometers using the data collected during the first step and the second step.

14. The method of claim 11, wherein collecting data from the at least two magnetometers of the inertial measurement unit while the inertial measurement unit is rigidly attached to the ferromagnetic object and repositioned includes repositioning the ferromagnetic object for at least one revolution around each Cartesian axis.

15. The method of claim 11, wherein calibrating the at least two magnetometers to account for ferromagnetic distortion includes centering the data collected from the at least two magnetometers.

16. The method of claim 15, wherein calibrating the at least two magnetometers to account for ferromagnetic distortion includes at least one of using a point correspondence transformation equation and using a geometric scaling and compensation equation.

17. The method of claim 11, wherein calibrating the at least two magnetometers to account for ferromagnetic distortion includes applying adjusted hard and soft iron compensation parameters.