Total Least Squares In-Field Identification for MEMS-Based Inertial Measurement Units ^†

Abstract

Inertial Measurement Units are widely used in various applications and, hardware-wise, they primarily consist of a tri-axial accelerometer and a tri-axial gyroscope. For low-end commercial employments, the low cost of the device is crucial: this makes MEMS-based sensors a popular choice in this context. However, MEMS-based transducers are prone to significant, non-uniform and environmental-condition-dependent systematic errors, that require frequent re-calibration to be eliminated. To this end, identification methods that can be performed in-field by non-expert users, without the need for high-precision or costly equipment, are of particular interest. In this paper, we propose an in-field identification procedure based on the Total Least Squares method for both tri-axial accelerometers and gyroscopes. The proposed identification model is linear and requires no prior knowledge of the parameters to be identified. It enables accelerometer calibration without the need for specific reference surface orientation relative to Earth’s gravity and allows gyroscope calibration to be performed independently of accelerometer data, without requiring the sensor’s sensitive axes to be aligned with the rotation axes during calibration. Experiments conducted on NXP sensors FXOS8700CQ and FXAS21002 demonstrated that using parameters identified by our method reduced cross-validation standard deviations by about two orders of magnitude compared to those obtained using manufacturer-provided parameters. This result indicates that our method enables the effective calibration of IMU sensor parameters, relying only on simple 3D-printed equipment and significantly improving IMU performance at minimal cost.

1. Introduction

Inertial Measurement Units (IMUs) are widely used in many applications where the position and attitude of a body need to be evaluated: civilian navigation systems [1], robotics [2], aerospace [3], medical applications [4], as well as commercial low-end portable devices such as computers or smartphones [5]. Six-axis IMUs are essentially composed of an accelerometer and a gyroscope: these are transducers that output analog or digital signals based on physical acceleration or angular velocity input, respectively. In both cases, a linear function can be used to describe the relationship between input and output, allowing us to determine a mathematical model for the sensor [6].

To meet the widespread demand for consumer-grade devices, IMUs based on Micro Electro-Mechanical Systems (MEMS-based) are of particular interest due to their advantages in cost, size, weight, and power consumption compared to mechanical or optical alternatives [7]. However, MEMS manufacturing technologies produce sensors having non-uniform model parameters, which the manufacturer only estimates on batches [6]. For these two reasons, if manufacturer-provided parameters are used in the sensors’ mathematical model, slight to large measurement errors can be expected. Additionally, if those measurements need to be post-processed (as for position and attitude estimation) the errors are likely to be further propagated and amplified [8].

Sensor errors are typically categorized into systematic and random errors. While random errors can only be accounted for using stochastic models [9], systematic errors are deterministic, i.e., they can be minimized using a calibration procedure: essentially, sensor model parameters are identified by comparing their mathematical model output with a known reference. Calibration techniques can be classified into two main types: Laboratory and in-field ones. While the former requires using high-precision and costly equipment, such as mechanical rotating tables, to provide a known reference, the latter can be conducted by the end user without the need for high-end apparatus or external sensors [6]. The low-repeatability of the manufacturing process, the low cost and the environmental-condition-dependent behavior that characterize MEMS-based sensors make in-field techniques a sensible choice for their calibration.

Accelerometer in-field calibration methods leverage the fact that, when the sensor is kept static, the Euclidean norm of the acceleration measured by the sensor in its reference frame (${\mathbf{a}}_m$) should be constant and equal to the one of the gravitational acceleration ($\mathbf{g}$). Most of these methods involve re-orienting the device N times to collect measurements ${\mathbf{a}}_{m,i}$, where $i=1,\cdots ,N$: the array collecting the sensor model parameters (${\mathit{\theta }}_a$) is then obtained by minimizing the cost function

$$C_a\left({\mathit{\theta }}_a\right)=\sum _{i=1}^N(\parallel {\mathbf{a}}_{m,i}\left({\mathit{\theta }}_a\right){\parallel }^2-{\parallel \mathbf{g}\parallel }^2{)}^2,$$

(1)

where $\parallel ·\parallel$ denotes the Euclidean norm of a vector. The problem is usually solved using Non-Linear Least Squares (NLS) methods [10,11,12]. Regarding gyroscopes, the Earth’s rotation rate cannot serve as a reliable reference due to its weak magnitude, which would result in an excessively low Signal to Noise Ratio (SNR) [5,6]. For this reason, accelerometer-dependent gyroscope calibration is often employed [11,12]. In this approach, the sensor is positioned in N static orientations by performing $N-1$ rotations: when the j-th non-static period begins, the gravity vector measured by the accelerometer (${\mathbf{g}}_{0,j}$) is fed to an algorithm that, by integrating gyroscope measurements acquired during sensor rotation, gives the output as an estimation of the gravity vector at the end of the rotation (${\mathbf{g}}_{g,j}$) that depends on the gyroscope’s acquired data. The array ${\mathit{\theta }}_g$ collecting the gyroscope model parameters is then obtained by minimizing the cost function

$$C_g\left({\mathit{\theta }}_g\right)=\sum _{j=1}^{N-1}{\parallel {\mathbf{g}}_{a,j}-{\mathbf{g}}_{g,j}\left({\mathit{\theta }}_g\right)\parallel }^2,$$

(2)

where ${\mathbf{g}}_{a,j}$ is the accelerometer-measured gravity vector at the end of the j-th rotation. These methods do not require placing the sensor in specific orientations, which enhances usability, but they suffer from two main drawbacks: (i) they rely on NLS optimization algorithms and consequently require a tentative initial parameter array, and (ii) the calibration outcome could be poor depending on the random orientations selected by the user. Moreover, if accelerometer-dependent gyroscope calibration is performed, any possible accelerometer calibration error is likely to be further propagated and amplified.

Methods for linear and non-iterative in-field calibration, which do not require any previous knowledge of the parameters to be calibrated, can be classified as identification methods and are highly desirable. For accelerometers, such methods were proposed in [13,14], but they rely on strict assumptions, including (i) the availability of a surface orthogonal to Earth’s gravity, and (ii) the sensor’s parameters being within known bounds. Other types of identification methods [15,16,17] for these devices are based on a suitable non-linear change of variables in the sensor model, which reformulates the classical NLS methods of [10,11,12] into linear ones. On the other hand, once the coefficients are determined, the inverse transformation of the non-linear change of variables needs to be performed, numerically propagating the errors in the estimated parameters. For gyroscopes, non-iterative methods have been proposed in [10,18]. However, these methods require rotating the sensor around its sensitive axes and relying on a mean angular velocity approach to compute rotations, rather than properly integrating gyroscope measurements. Assuming that the rotation axis coincides with one of the sensor’s axes and using the mean angular velocity for integration are weak assumptions, especially given that the calibration is intended to be performed in-field by the end user. In fact, if the rotations are to be performed manually, (i) precisely rotating around one of the sensor’s sensitive axes can be extremely challenging and (ii) achieving a smooth angular velocity profile can be difficult. Therefore, assuming an unknown rotation axis and performing a small time-step integration of gyroscope measurements would be desirable approaches to improve the robustness of the calibration. Moreover, these methods rely on accelerometer measurements to estimate reference rotation angles, which can result in error propagation. Additionally, they do not account for cross-sensitivities between the sensor axes.

In this paper, we propose an in-field identification method that extends the approach presented in [19] for accelerometers to include gyroscope calibration as well. This extended method is based on (i) a predetermined sequence of rotations used to acquire sensor outputs, and (ii) the solution of two Total Least Squares (TLS) problems. The sensor is placed in a 3D-printed multifaceted housing, which is rotated with respect to a sloped surface used as a reference. Each rotation is manually performed around a fixed axis that, in general, is assumed not to coincide with any of the sensor triad axes. Small time-step trapezoidal integration of gyroscope measurements is employed to obtain a reliable estimation of rotation angles. These estimations are then compared with the rotation angles expected given the geometry of the IMU’s housing and the chosen rotation sequence, allowing for accelerometer-independent gyroscope calibration. The TLS approach ensures that potential errors in the rotations are accounted for during calibration. In this manner, sensor model parameters, including sensitivities, cross-sensitivities and bias terms, are identified. In summary, our method’s contributions are as follows:

• Our method does not require the calibration reference to have a known orientation with respect to Earth’s gravity, addressing a key shortcoming of existing accelerometer identification methods [13,14]. This enhances the in-field usability of our method.
• In the gyroscope calibration model, the axis of each rotation is assumed to be unknown, i.e., non-coincident with any of the sensor’s sensitive axes. This contrasts with state-of-the-art gyroscope identification methods [10,18] and improves the robustness of the identification, especially since rotations are performed manually.
• Accelerometer and gyroscope calibration are independent of each other, meaning that any errors in accelerometer calibration do not propagate to gyroscope calibration. This solves one of the drawbacks of [10,11,12,18].
• Both accelerometer and gyroscope identification models are linear, eliminating the need for any initial guess regarding sensor parameters. In NLS-based methods [11,12], a poor choice of the initial guess can lead to optimization being trapped in a local minimum of the cost function, resulting in sub-optimal calibration outcomes.
• The TLS approach is used to solve the linear identification systems of equations, allowing for the consideration of potential errors in the data acquisition process. In contrast, other linear methods like [10,18] use least squares solutions that assume the measurements acquired during the data collection process to be error-free, which is often not the case when rotations for in-field calibration are performed manually.

The rest of the paper is organized as follows: Section 2 presents the description of the identification algorithms, while experiments performed using NXP sensors FXOS8700CQ and FXAS21002 mounted on a single chip are presented and discussed in Section 3. Conclusions are drawn in Section 4.

2. Identification Algorithms Description

The proposed identification method requires placing a “single-board” IMU—namely, a device where the accelerometer and the gyroscope are rigidly mounted on the same electronic board and mechanical frame—inside a multifaceted 3D-printed housing, and re-orient such housing into N different orientations with respect to (w.r.t.) a reference surface. In practice, this is accomplished by positioning the housing on each of its facets according to a predetermined sequence of $N-1$ rotations. Let $Oxyz$ be a fixed reference frame rigidly attached to the sloped surface where the calibration takes place and $O^{\prime }{x}^{\prime }{y}^{\prime }{z}^{\prime }$ be a mobile reference frame rigidly attached to the IMU.

We will refer to them as the fixed frame (F-frame) and the mobile frame (M-frame), respectively, and we assume accelerometer and gyroscope orthogonal triads to be coincident with the M-frame. When the calibration procedure starts, the M-frame is assumed to be randomly tilted w.r.t. the F-frame (Figure 1), meaning that the IMU is placed in an unknown orientation w.r.t. its housing. In the following subsections, accelerometer and gyroscope sensor models are introduced and their calibration algorithms are described.

2.1. Accelerometer

The mathematical model for a tri-axial accelerometer, which correlates the acceleration perceived by the sensor expressed w.r.t. the M-frame ${\left\{\mathbf{a}\right\}}_M\in {\mathbb{R}}_{3\times 1}$ with the sensor’s raw output $\mathbf{v}={[v_1,v_2,v_3]}^T$ is given by [6]:

$$\mathbf{v}={\mathbf{A}}^{\prime }{\left\{\mathbf{a}\right\}}_M+{\mathbf{b}}^{\prime }+{\mathit{\zeta }}_a^{\prime },$$

(3)

where ${\mathbf{A}}^{\prime }\in {\mathbb{R}}_{3\times 3}$ is a scale factor matrix, ${\mathbf{b}}^{\prime }\in {\mathbb{R}}_{3\times 1}$ is an array collecting bias terms and ${\mathit{\zeta }}_a^{\prime }\in {\mathbb{R}}_{3\times 1}$ contains measurement noise. Although this model accurately represents the physical behavior of the sensor, its inverse form is preferable for practical uses:

$${\left\{\mathbf{a}\right\}}_M=\mathbf{A}(\mathbf{v}+{\mathit{\zeta }}_a)+\mathbf{b},$$

(4)

where $\mathbf{A}={\mathbf{A}\prime }^{-1},\mathbf{b}=-{\mathbf{A}\prime }^{-1}{\mathbf{b}}^{\prime }$ and ${\mathit{\zeta }}_a=-{\mathit{\zeta }}_a^{\prime }$. The identification procedure aims to determine the elements of $\mathbf{A}$ and $\mathbf{b}$, so that raw sensor output $\mathbf{v}$ can be transformed into an acceleration measurement ${\left\{\mathbf{a}\right\}}_M$. $\mathbf{A}$ is assumed to be symmetric and noise terms to be white Gaussian with zero mean [6]. We call diagonal elements $A_{pp}$ of matrix $\mathbf{A}$sensitivities, while its off-diagonal elements $A_{pq}$ are referred to as cross-sensitivities. Cross-sensitivities are non-zero when non-orthogonality in the sensor’s triad arises due to manufacturing tolerances, and they also account for cross-talk effect between different channels caused by the sensor’s electronics. These non-idealities are particularly relevant for MEMS-based sensors [20] and our identification model takes them into account. We assume the gravitational acceleration to have a constant magnitude equal to g that we denote as ${\{·\}}_{M,i}$; the expression of a vector in M-frame coordinates when the IMU is placed in the i-th orientation of the procedure, for $i=1,\cdots ,N$. Accordingly, the acceleration perceived by the sensor in the first orientation of the procedure can be written as follows:

$${\left\{\mathbf{a}\right\}}_{M,1}=g{\left\{\mathbf{n}\right\}}_{M,1},$$

(5)

where ${\left\{\mathbf{n}\right\}}_{M,1}={[n_1,n_2,n_3]}^T$ is a unit vector (i.e., its Euclidean norm is one) expressing Earth’s gravity direction in the mobile frame. ${\left\{\mathbf{n}\right\}}_{M,1}$ is generally unknown, given that the reference is randomly oriented w.r.t. Earth’s gravity, and it will be determined as a sub-product of the identification. If the IMU is re-oriented in N different orientations during the calibration procedure, the acceleration perceived by the sensor in the i-th orientation for $i=1,\cdots ,N$ can be expressed as follows:

$${\left\{\mathbf{a}\right\}}_{M,i}={\mathbf{R}}_i{\left\{\mathbf{a}\right\}}_{M,1},$$

(6)

where ${\mathbf{R}}_i$ is a rotation matrix defined by the known sequence of rotations performed w.r.t. the F-frame. If we substitute Equation (4) into Equation (6), we can correlate the raw output of the sensor in the i-th configuration with the acceleration in the initial configuration as follows:

$${\left\{\mathbf{a}\right\}}_{M,i}=\mathbf{A}({\mathbf{v}}_i+{\mathit{\zeta }}_{a,i})+\mathbf{b}={\mathbf{R}}_i{\left\{\mathbf{a}\right\}}_{M,1}.$$

(7)

Since measurement noise can be well-modeled as white Gaussian with zero mean, the average of k raw output samples acquired in the same orientation is considered in order to mitigate noise effects, namely:

$${\overline{\mathbf{v}}}_i={\displaystyle \frac{{\sum }_{s=1}^k({\mathbf{v}}_{i,s}+{\mathit{\zeta }}_{a,i,s})}{k}},k\gg 1.$$

(8)

In this way, the noiseless mathematical model is obtained as follows:

$${\left\{\mathbf{a}\right\}}_{M,i}=\mathbf{A}{\overline{\mathbf{v}}}_i+\mathbf{b}={\mathbf{R}}_i{\left\{\mathbf{a}\right\}}_{M,1}.$$

(9)

Re-arranging terms of Equation (9) so as to highlight known quantities and parameters to be calibrated yields the following:

$${\mathbf{M}}_i\mathit{\theta }={\mathbf{0}}_{3\times 1},$$

(10)

where

$${\mathbf{M}}_i=\left[\begin{array}{cccccc}{\overline{v}}_{i,1}& {\overline{v}}_{i,2}& {\overline{v}}_{i,3}& 0& 0& 0\\ 0& {\overline{v}}_{i,1}& 0& {\overline{v}}_{i,2}& {\overline{v}}_{i,3}& 0& {\mathbf{I}}_{3\times 3}& -{\mathbf{R}}_i\\ 0& 0& {\overline{v}}_{i,1}& 0& {\overline{v}}_{i,2}& {\overline{v}}_{i,3}\end{array}\right]$$

(11)

and

$$\mathit{\theta }={[A_{11},A_{12},A_{13},A_{22},A_{23},A_{33},b_1,b_2,b_3,n_1,n_2,n_3]}^T.$$

(12)

We denoted as ${\mathbf{I}}_{3\times 3}$ the third order identity matrix, as $A_{uw}$ the element in the u-th row and w-th column of matrix $\mathbf{A}$ and as $b_1,b_2,b_3$ the elements of vector $\mathbf{b}$. According to our identification method, the gravity direction in the start configuration ${\left\{\mathbf{n}\right\}}_{M,1}$, expressed in the mobile frame, is also unknown and included in the parameters to be determined. By stacking Equation (10) for $i=1,\cdots ,N$ measurements, we obtain a homogeneous system of $3N$ linear equations in 12 unknowns, which is over-determined if $N>4$:

$$\left[\begin{array}{c}{\mathbf{M}}_1\\ ⋮\\ {\mathbf{M}}_N\end{array}\right]\mathit{\theta }=\mathbf{M}\mathit{\theta }={\mathbf{0}}_{3N\times 1},\mathbf{M}\in {\mathbb{R}}_{3N\times 12}.$$

(13)

Equation (13) can be solved for $\mathit{\theta }$ with the linear TLS method (see [21] for the theoretical foundations, or [22,23] for practical applications). Using the TLS method allows us to account for possible errors in the elements of matrix $\mathbf{M}$ originated in the data acquisition procedure. Please note that $\mathit{\theta }$ contains the actual parameters needed by the model, and not a combination of them, which is the second shortcoming of previous works [15,16,17]. In ideal conditions—namely, without any measurement noise or model error—$\mathbf{M}$ should be full-column rank deficient—namely, $rank\left(\mathbf{M}\right)<12$—and a solution for Equation (13) different from the trivial one (which is $\mathit{\theta }=\mathbf{0}$) could be found. In real applications, it is necessary to address a compatible system that is closest to Equation (13) in terms of the Frobenious norm:

$$\widehat{\mathbf{M}}\widehat{\mathit{\theta }}=\mathbf{0},$$

(14)

where $\widehat{\mathit{\theta }}$, which is the solution to Equation (13), contains the unknown parameters estimations and $\widehat{\mathbf{M}}$ has rank 11 and minimizes $\parallel \widehat{\mathbf{M}}{-\mathbf{M}\parallel }_F$, where ${\parallel ·\parallel }_F$ is the Frobenius norm. Matrix $\widehat{\mathbf{M}}$ can be computed by means of the "economy size" singular value decomposition of $\mathbf{M}$, given by the following:

$$\mathbf{M}=\mathbf{U}\left[\begin{array}{c}diag\left({\lambda }_l\right)\\ \mathbf{0}\end{array}\right]{\mathbf{V}}^T,$$

(15)

where $\mathbf{U}\in {\mathbb{R}}_{3N\times 3N}$ and $\mathbf{V}\in {\mathbb{R}}_{12\times 12}$ are orthonormal matrices and $diag\left({\lambda }_l\right)\in {\mathbb{R}}_{12\times 12}$ is a diagonal matrix containing the singular values ${\lambda }_l$, for $l=1,\cdots ,12$ of $\mathbf{M}$ sorted in decreasing order. $\widehat{\mathbf{M}}$ can be obtained as follows:

$$\widehat{\mathbf{M}}=\mathbf{M}-{\lambda }_{12}{\mathbf{U}}_{12}{\mathbf{V}}_{12}^T,$$

(16)

being ${\lambda }_{12}$ the smallest singular value of $\mathbf{M}$ and ${\mathbf{U}}_{12},{\mathbf{V}}_{12}^T$ the 12-th column of $\mathbf{U}$ and $\mathbf{V}$, respectively. A solution to Equation (14) is given by $\widehat{\mathit{\theta }}={\mathbf{V}}_{12}$. The solution to the proposed identification problem is obtained by normalizing $\widehat{\mathit{\theta }}$ so that its The last three elements have unitary Euclidean norm. This generates an ambiguity regarding the sign of the solution. However, the ambiguity can be removed by making an ex-post check on the sign of the last element of $\mathit{\theta }$, for example, if in the initial orientation of the procedure $z^{\prime }$ is roughly pointing in the opposite direction w.r.t.$\mathbf{g}$, then the last element of $\mathit{\theta }$ will be imposed to have negative sign), namely:

$$\mathit{\theta }={\displaystyle \frac{\widehat{\mathit{\theta }}}{\parallel {\widehat{\mathit{\theta }}}_{10:12}\parallel }},$$

(17)

where ${\widehat{\mathit{\theta }}}_{10:12}$ is the array collecting the last three elements of $\widehat{\mathit{\theta }}$.

Percentage relative standard deviation estimations for the identified parameters can be used to infer the accuracy of the obtained results. The estimation is performed according to the TLS techniques, and by assuming $n_3$ as the element having the maximum accuracy. Accordingly, the covariance matrix of the TLS solution may be approximated as [21]:

$${\mathbf{C}}_{\mathit{\theta }}={\widehat{\sigma }}_M^2(1+\parallel {\mathit{\theta }}_{1:11}{\parallel }^2){\left({\widehat{\mathbf{M}}}_{1:11}^T{\widehat{\mathbf{M}}}_{1:11}\right)}^{-1},$$

(18)

where

$${\sigma }_M={\displaystyle \frac{{\lambda }_{12}}{\sqrt{3N-12}}}$$

(19)

and ${\widehat{\mathbf{M}}}_{1:11}$ is obtained by discarding the last column of $\widehat{\mathbf{M}}$, whereas ${\mathit{\theta }}_{1:11}$ is the array containing the first 11 elements of $\mathit{\theta }$. We finally obtain the percentage relative standard deviation estimation as [23]:

$${\sigma }_{\%,{\mathit{\theta }}_l}=100{\displaystyle \frac{\sqrt{{\mathbf{C}}_{\mathit{\theta }}(l,l)}}{|{\mathit{\theta }}_l|}},l=1,\cdots ,11;$$

(20)

where ${\mathbf{C}}_{\mathit{\theta }}(l,l)$ is the l-th diagonal term of ${\mathbf{C}}_{\mathit{\theta }}$ and ${\mathit{\theta }}_l$ the l-th element of $\mathit{\theta }$.

2.2. Gyroscope

Analogous to Section 2.1, we consider the angular velocity perceived by the gyroscope ${\left\{\omega \right\}}_M$, expressed in the M-frame, to be related to the sensor raw output $\mathbf{r}={[r_1,r_2,r_3]}^T$ by the equation:

$${\left\{\omega \right\}}_M=\mathbf{G}\mathbf{r}+\mathbf{d}+{\mathit{\zeta }}_g,$$

(21)

where diagonal elements $G_{pp}$ of the symmetric matrix $\mathbf{G}\in {\mathbb{R}}_{3\times 3}$ are referred to as sensitivities, off-diagonal ones $G_{pq}$ are cross-sensitivities (which are, in general, non-zero), $\mathbf{d}={[d_1,d_2,d_3]}^T$ collects bias terms and ${\mathit{\zeta }}_g\in {\mathbb{R}}_{3\times 1}$ contains measurement noise. Previous identification methods [10,18] consider the cross-sensitivities to be negligible and this can lead to inaccurate calibration for MEMS-based sensors [20].

Let ${\mathbf{u}}_x,{\mathbf{u}}_y,{\mathbf{u}}_z$ be three unit vectors that represent F-frame axes directions (Figure 1). When the sensor is placed in the first orientation of the calibration procedure, by using the Euler–Rodrigues formula [24], we can express ${\mathbf{u}}_x,{\mathbf{u}}_y,{\mathbf{u}}_z$ in M-frame coordinates as follows:

$$\begin{array}{c}\hfill {\left\{{\mathbf{u}}_x\right\}}_{M,1}=cos\varphi \left[\begin{array}{c}1\\ 0\\ 0\end{array}\right]+sin\varphi (\mathbf{e}\times \left[\begin{array}{c}1\\ 0\\ 0\end{array}\right])+(1-cos\varphi )\mathbf{e}(\mathbf{e}·\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right]),\\ \hfill {\left\{{\mathbf{u}}_y\right\}}_{M,1}=cos\varphi \left[\begin{array}{c}0\\ 1\\ 0\end{array}\right]+sin\varphi (\mathbf{e}\times \left[\begin{array}{c}0\\ 1\\ 0\end{array}\right])+(1-cos\varphi )\mathbf{e}(\mathbf{e}·\left[\begin{array}{c}0\\ 1\\ 0\end{array}\right]),\\ \hfill {\left\{{\mathbf{u}}_z\right\}}_{M,1}=cos\varphi \left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]+sin\varphi (\mathbf{e}\times \left[\begin{array}{c}0\\ 0\\ 1\end{array}\right])+(1-cos\varphi )\mathbf{e}(\mathbf{e}·\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]),\end{array}$$

(22)

where $\mathbf{e}={[e_1,e_2,e_3]}^T$ and $\varphi$ are the axis and the angle that define the rotation needed to transform the M-frame into the F-frame, while $·,\times$ denote the Euclidean dot and cross products, respectively. In this first IMU orientation, we assume misalignments between the M-frame and the F-frame to be small (i.e., $\varphi \approx 0$), so that Equation (22) can be approximated as:

$$\begin{array}{c}\hfill {\left\{{\mathbf{u}}_x\right\}}_{M,1}\approx \left[\begin{array}{c}1\\ \varphi e_3\\ -\varphi e_2\end{array}\right],\\ \hfill {\left\{{\mathbf{u}}_y\right\}}_{M,1}\approx \left[\begin{array}{c}-\varphi e_3\\ 1\\ \varphi e_1\end{array}\right],\\ \hfill {\left\{{\mathbf{u}}_z\right\}}_{M,1}\approx \left[\begin{array}{c}\varphi e_2\\ -\varphi e_1\\ 1\end{array}\right].\end{array}$$

(23)

According to our calibration method, the IMU needs to be rotated w.r.t. F-frame fixed axes; for each j-th rotation, with $j=1,\cdots ,N-1$, we can express in the M-frame the angular velocity ${\left\{\omega \left(t\right)\right\}}_{M,j}$ perceived by the gyroscope at time instant t as follows:

$${\left\{\omega \left(t\right)\right\}}_{M,j}=\dot{\delta }\left(t\right){\left\{{\mathbf{u}}_h\right\}}_{M,j},$$

(24)

where $\dot{\delta }\left(t\right)$ is the angular velocity magnitude at time instant t and $h=x,y,z$ depending on whether the rotation is performed around $x,y$ or z axis, respectively. Please notice that ${\mathbf{u}}_x,{\mathbf{u}}_y,{\mathbf{u}}_z$ are fixed, namely they remain the same throughout the calibration procedure, but, in general, their expressions in the gyroscope triad (the M-frame) change between the different orientations of the procedure. We can relate the expression of ${\mathbf{u}}_h$ in the M-frame to ${\mathbf{u}}_h$ assumed in the initial configuration as follows:

$${\left\{{\mathbf{u}}_h\right\}}_{M,j}={\mathbf{R}}_j{\left\{{\mathbf{u}}_h\right\}}_{M,1},$$

(25)

where ${\mathbf{R}}_j$ is the rotation matrix that transforms a vector expressed in M,1-frame coordinates into a vector expressed in M,j-frame ones. Such a rotation matrix is coincident with the one introduced in Equation (6) and its results are determined by the known sequence of rotations performed. By substituting Equations (25) and (21) in Equation (24), we obtain:

$${\left\{\omega \left(t\right)\right\}}_{M,j}=\mathbf{G}\mathbf{r}\left(t\right)+\mathbf{d}+{\mathit{\zeta }}_g\left(t\right)=\dot{\delta }\left(t\right){\mathbf{R}}_j{\left\{{\mathbf{u}}_h\right\}}_{M,1}.$$

(26)

Measurement noise can be removed by using proper filtering techniques: analogic or digital filters can be embedded directly in the IMU electronic board or implemented in the data acquisition system if the sampling time is sufficiently small. The noiseless calibration model can then be written as follows:

$${\left\{\omega \left(t\right)\right\}}_{M,j}=\mathbf{G}\mathbf{r}\left(t\right)+\mathbf{d}=\dot{\delta }\left(t\right){\mathbf{R}}_j{\left\{{\mathbf{u}}_h\right\}}_{M,1}.$$

(27)

By integrating the right-most term of Equation (27) in the time interval $[t_{0,j},t_{1,j}]$ where the j-th rotation takes place, we obtain:

$${\int }_{t_{0,j}}^{t_{1,j}}\dot{\delta }\left(t\right){\mathbf{R}}_j{\left\{{\mathbf{u}}_h\right\}}_{M,1}dt={\mathrm{\Delta }}_j{\mathbf{R}}_j{\left\{{\mathbf{u}}_h\right\}}_{M,1},$$

(28)

where ${\mathrm{\Delta }}_j$ is the j-th expected rotation angle, whose numerical value depends on the geometry of IMU’s 3D-printed housing and the chosen rotation sequence. Since the estimation of ${\mathrm{\Delta }}_j$ does not rely on accelerometer measurements (as was the case in [10,18]), the results of gyroscope calibration are independent of those of the accelerometer. Integrating in the same time interval the central term of Equation (27) one obtains:

$${\int }_{t_{0,j}}^{t_{1,j}}(\mathbf{G}\mathbf{r}\left(t\right)+\mathbf{d})dt=\left[\begin{array}{c}{G}_{11}{H}_{1,j}+G_{12}{H}_{2,j}+G_{13}{H}_{3,j}+T_j{d}_1\\ G_{12}{H}_{1,j}+G_{22}{H}_{2,j}+G_{23}{H}_{3,j}+T_j{d}_2\\ G_{13}{H}_{1,j}+G_{23}{H}_{2,j}+G_{33}{H}_{3,j}+T_j{d}_3\end{array}\right],$$

(29)

where

$$\begin{array}{c}\hfill H_{1,j}={\int }_{t_{0,j}}^{t_{1,j}}{r}_1\left(t\right)dt,\\ \hfill H_{2,j}={\int }_{t_{0,j}}^{t_{1,j}}{r}_2\left(t\right)dt,\\ \hfill H_{3,j}={\int }_{t_{0,j}}^{t_{1,j}}{r}_3\left(t\right)dt,\\ \hfill T_j=t_{1,j}-t_{0,j},\end{array}$$

$G_{uw}$ is the element in the u-th row and w-th column of matrix $\mathbf{G}$ and $d_1,d_2,d_3$ the elements of vector $\mathbf{d}$. As introduced in Section 1, $H_{1,j},H_{2,j},H_{3,j}$ are obtained through the integration of gyroscope raw measurements and not by using a mean value approach. Combining Equations (27)–(29) we obtain:

$$\left[\begin{array}{c}{G}_{11}{H}_{1,j}+G_{12}{H}_{2,j}+G_{13}{H}_{3,j}+T_j{d}_1\\ G_{12}{H}_{1,j}+G_{22}{H}_{2,j}+G_{23}{H}_{3,j}+T_j{d}_2\\ G_{13}{H}_{1,j}+G_{23}{H}_{2,j}+G_{33}{H}_{3,j}+T_j{d}_3\end{array}\right]-{\mathrm{\Delta }}_j{\mathbf{R}}_j{\left\{{\mathbf{u}}_h\right\}}_{M,1}={\mathbf{0}}_{3\times 1}.$$

(30)

Please notice that Equation (30) can assume a different form depending on weather the j-th rotation is performed around $x,y$ or z, namely:

$$\begin{array}{c}\hfill \left[\begin{array}{c}{G}_{11}{H}_{1,j}+G_{12}{H}_{2,j}+G_{13}{H}_{3,j}+T_j{d}_1\\ G_{12}{H}_{1,j}+G_{22}{H}_{2,j}+G_{23}{H}_{3,j}+T_j{d}_2\\ G_{13}{H}_{1,j}+G_{23}{H}_{2,j}+G_{33}{H}_{3,j}+T_j{d}_3\end{array}\right]-{\mathrm{\Delta }}_j{\mathbf{R}}_j{\left\{{\mathbf{u}}_x\right\}}_{M,1}={\mathbf{0}}_{3\times 1},\mathrm{if}h=x\\ \hfill \left[\begin{array}{c}{G}_{11}{H}_{1,j}+G_{12}{H}_{2,j}+G_{13}{H}_{3,j}+T_j{d}_1\\ G_{12}{H}_{1,j}+G_{22}{H}_{2,j}+G_{23}{H}_{3,j}+T_j{d}_2\\ G_{13}{H}_{1,j}+G_{23}{H}_{2,j}+G_{33}{H}_{3,j}+T_j{d}_3\end{array}\right]-{\mathrm{\Delta }}_j{\mathbf{R}}_j{\left\{{\mathbf{u}}_y\right\}}_{M,1}={\mathbf{0}}_{3\times 1},\mathrm{if}h=y\\ \hfill \left[\begin{array}{c}{G}_{11}{H}_{1,j}+G_{12}{H}_{2,j}+G_{13}{H}_{3,j}+T_j{d}_1\\ G_{12}{H}_{1,j}+G_{22}{H}_{2,j}+G_{23}{H}_{3,j}+T_j{d}_2\\ G_{13}{H}_{1,j}+G_{23}{H}_{2,j}+G_{33}{H}_{3,j}+T_j{d}_3\end{array}\right]-{\mathrm{\Delta }}_j{\mathbf{R}}_j{\left\{{\mathbf{u}}_z\right\}}_{M,1}={\mathbf{0}}_{3\times 1},\mathrm{if}h=z\end{array}$$

(31)

A relevant contribution to gyroscope identification methodology introduced in this work can now be highlighted: the simplification introduced by Equation (23) allows the inclusion of the expressions for the rotation axes in the M-frame at the beginning of the procedure (namely ${\left\{{\mathbf{u}}_x\right\}}_{M,1},{\left\{{\mathbf{u}}_y\right\}}_{M,1},{\left\{{\mathbf{u}}_z\right\}}_{M,1}$) into the identification model without producing any non-linearity into the equations and without assuming that ${\left\{{\mathbf{u}}_x\right\}}_{M,1}={[1,0,0]}^T,$${\left\{{\mathbf{u}}_y\right\}}_{M,1}={[0,1,0]}^T$ or ${\left\{{\mathbf{u}}_z\right\}}_{M,1}={[0,0,1]}^T$ (i.e., without assuming alignment of the rotation axes with the gyroscope’s sensitive axes). In this way, one of the shortcomings of previous gyroscope identification methods [10,18] is solved. One may argue that the simplification introduced in Equation (23) is valid only if $\varphi$ (which is unknown and included, together with $\mathbf{e}$, in the parameters to be determined) is sufficiently small: however, this assumption is justified for industrially manufactured sensors, and is further confirmed by the experimental results shown in Section 3. By substituting the expressions of ${\left\{{\mathbf{u}}_x\right\}}_{M,1},{\left\{{\mathbf{u}}_y\right\}}_{M,1},{\left\{{\mathbf{u}}_z\right\}}_{M,1}$ with those given by Equation (23) and re-organizing terms of Equation (31) so as to highlight known quantities and unknown parameters to be identified, each j-th rotation can be associated with a linear system of equations of the form:

$${\mathbf{L}}_j\mathit{\vartheta }={\mathbf{0}}_{3\times 1},$$

(32)

where

$$\mathit{\vartheta }={[G_{11},G_{12},G_{13},G_{22},G_{23},G_{33},d_1,d_2,d_3,\varphi e_3,\varphi e_2,-\varphi e_1,1]}^T$$

(33)

is the gyroscope-related counterpart of $\mathit{\theta }$, introduced in Equation (12) to collect accelerometer-related parameters to be identified: analogously, $\mathit{\vartheta }$ contains gyroscope model parameters, quantities related to $\mathbf{e},\varphi$ and a known term, i.e., its last element ${\mathit{\vartheta }}_{13}=1$. ${\mathbf{L}}_j$ can be expressed as composition of two different matrices ${\mathbf{L}}_{j,1}$ and ${\mathrm{\Delta }}_j{\mathbf{L}}_{j,2}$, namely:

$${\mathbf{L}}_j=\left[\begin{array}{cc}{\mathbf{L}}_{j,1}& {\mathrm{\Delta }}_j{\mathbf{L}}_{j,2}\end{array}\right],$$

(34)

where ${\mathbf{L}}_{j,1}$ is:

$${\mathbf{L}}_{j,1}=\left[\begin{array}{cccccc}{H}_{1,j}& H_{2,j}& H_{3,j}& 0& 0& 0\\ 0& H_{1,j}& 0& H_{2,j}& H_{3,j}& 0& T_j{\mathbf{I}}_{3\times 3}\\ 0& 0& H_{1,j}& 0& H_{2,j}& H_{3,j}\end{array}\right]$$

(35)

and has the same form independent of the axis around which the j-th rotation is performed, while ${\mathbf{L}}_{j,2}$ is different depending on weather the j-th rotation is performed around $x,y$ or z axis and can be defined as follows:

$${\mathbf{L}}_{j,2}=\left\{\begin{array}{cc}\left[\begin{array}{cccc}\mp {\mathbf{R}}_{j,2}& \pm {\mathbf{R}}_{j,3}& {\mathbf{0}}_{3\times 1}& \mp {\mathbf{R}}_{j,1}\end{array}\right],\hfill & \mathrm{if}h=x\hfill \\ \left[\begin{array}{cccc}\pm {\mathbf{R}}_{j,1}& {\mathbf{0}}_{3\times 1}& \pm {\mathbf{R}}_{j,3}& ,\mp {\mathbf{R}}_{j,2}\end{array}\right],\hfill & \mathrm{if}h=y\hfill \\ \left[\begin{array}{cc}{\mathbf{0}}_{3\times 1}& \mp {\mathbf{R}}_j\end{array}\right],\hfill & \mathrm{if}h=z\hfill \end{array}\right.,$$

(36)

being ${\mathbf{R}}_{j,p}$ the p-th column of matrix ${\mathbf{R}}_j$, whose sign depends on weather the j-th rotation is performed counter-clockwise or clockwise, respectively. The mathematical derivation of Equation (32) can be found in Appendix A. By stacking Equation (32) for $j=1,\cdots ,N-1$, one obtains an homogeneous system of $3(N-1)$ linear equations in 13 unknowns, which is over-determined if $N\ge 6$:

$$\left[\begin{array}{c}{\mathbf{L}}_1\\ ⋮\\ {\mathbf{L}}_{N-1}\end{array}\right]\mathit{\vartheta }=\mathbf{L}\mathit{\vartheta }={\mathbf{0}}_{3(N-1)\times 1},\mathbf{L}\in {\mathbb{R}}_{3(N-1)\times 13}$$

(37)

Equation (37) can be solved for $\mathit{\vartheta }$ using the linear TLS method, following the same approach taken in Section 2.1 [21,22,23]. Even in this case, employing the TLS solution for Equation (37) ensures that potential errors in the elements of $\mathbf{L}$, arising from data acquisition inaccuracies, are accounted for. This represents another significant improvement over the methods proposed in [10,18], especially since manual rotations are likely to introduce inaccuracies in the measurement process. We address the compatible system closest to Equation (37) in terms of the Frobenius norm:

$$\widehat{\mathbf{L}}\widehat{\mathit{\vartheta }}=\mathbf{0},$$

(38)

where $\widehat{\mathit{\vartheta }}$ is the solution to Equation (37), containing estimated parameters, and $\widehat{\mathbf{L}}$ is such that $rank\left(\widehat{\mathbf{L}}\right)=12$ and $\parallel \widehat{\mathbf{L}}{-\mathbf{L}\parallel }_F$ is minimized. Exploiting the “economy size” singular value decomposition, $\mathbf{L}$ can be written as follows:

$$\mathbf{L}={\mathbf{U}}^{\prime }\left[\begin{array}{c}diag\left({\beta }_q\right)\\ \mathbf{0}\end{array}\right]{\mathbf{V}\prime }^T,$$

(39)

being ${\mathbf{U}}^{\prime }\in {\mathbb{R}}_{3(N-1)\times 3(N-1)}$ and ${\mathbf{V}}^{\prime }\in {\mathbb{R}}_{13\times 13}$ orthonormal matrices, while $diag\left({\beta }_q\right)\in {\mathbb{R}}_{13\times 13}$ is a diagonal matrix containing singular values ${\beta }_q$ of $\mathbf{L}$, for $q=1,\cdots ,13$, sorted in decreasing order. We can obtain $\widehat{\mathbf{L}}$ as follows:

$$\widehat{\mathbf{L}}=\mathbf{L}-{\beta }_{13}{\mathbf{U}}_{13}^{\prime }{\mathbf{V}}_{13}^{\prime T},$$

(40)

where ${\beta }_{13}$ is the smallest singular value of $\mathbf{L}$ and ${\mathbf{U}}_{13}^{\prime },{\mathbf{V}}_{13}^{\prime }$ are the 13-th column of ${\mathbf{U}}^{\prime }$ and ${\mathbf{V}}^{\prime }$, respectively. A solution to Equation (38) is $\widehat{\mathit{\vartheta }}={\mathbf{V}}_{13}^{\prime }$. The solution $\mathit{\vartheta }$ to the proposed identification problem can be obtained by taking into account that its last element ${\mathit{\vartheta }}_{13}$ must be equal to 1 according to Equation (33), namely:

$$\mathit{\vartheta }={\displaystyle \frac{\widehat{\mathit{\vartheta }}}{{\widehat{\mathit{\vartheta }}}_{13}}},$$

(41)

where ${\widehat{\mathit{\vartheta }}}_{13}$ is the last element of $\widehat{\mathit{\vartheta }}$. We can then obtain $\mathbf{e}$ and $\varphi$ as follows:

$$\mathbf{e}={\displaystyle \frac{{\mathit{\vartheta }}_{12:10}}{\parallel {\mathit{\vartheta }}_{12:10}\parallel }},\varphi =\parallel {\mathit{\vartheta }}_{12:10}\parallel ,$$

(42)

being ${\mathit{\vartheta }}_{12:10}={[-{\mathit{\vartheta }}_{12},{\mathit{\vartheta }}_{11},{\mathit{\vartheta }}_{10}]}^T$, where ${\mathit{\vartheta }}_q$ is the q-th element of $\mathit{\vartheta }$. Since ${\mathit{\vartheta }}_{13}=1$ is known, it is the element of $\mathit{\vartheta }$ having the maximum accuracy. Thus, we can approximate the covariance matrix of the TLS solution as follows:

$${\mathbf{C}}_{\mathit{\vartheta }}={\widehat{\sigma }}_L^2(1+\parallel {\mathit{\vartheta }}_{1:12}{\parallel }^2){\left({\widehat{\mathbf{L}}}_{1:12}^T{\widehat{\mathbf{L}}}_{1:12}\right)}^{-1},$$

(43)

where

$${\widehat{\sigma }}_L={\displaystyle \frac{{\beta }_{13}}{\sqrt{3(N-1)-13}}},$$

(44)

${\widehat{\mathbf{L}}}_{1:12}$ contains the first 12 columns of $\widehat{\mathbf{L}}$, while ${\mathit{\vartheta }}_{1:12}$ contains the first 12 elements of $\mathit{\vartheta }$. The percentage relative standard deviation is then obtained as follows:

$${\sigma }_{\%,{\mathit{\vartheta }}_q}=100{\displaystyle \frac{\sqrt{{\mathbf{C}}_{\mathit{\vartheta }}(q,q)}}{|{\mathit{\vartheta }}_q|}},q=1,\cdots ,12;$$

(45)

where ${\mathbf{C}}_{\mathit{\vartheta }}(q,q)$ is the q-th diagonal term of ${\mathbf{C}}_{\mathit{\vartheta }}$.

3. Experimental Procedure and Results

Simple 3D-printed IMU housings and references were produced to test the proposed method. This equipment includes a regular rectangular prism (prismatic) housing (Figure 2a), an 18-faced housing (Figure 2b) obtained by chamfering the edges of the former housing at 45°, a reference with three orthogonal surfaces and a base parallel to one of those surfaces (Figure 2c), and a reference with three orthogonal surfaces and a sloped base (Figure 2d). For each housing, a re-orientation scheme was defined to place the housing in at least every possible orientation where three of its faces are in contact with reference surfaces. Accordingly, $N=24$ and $N=74$ orientations were achieved in our tests for the prismatic and 18-faced housings, respectively.

Figure 3 illustrates a schematic of the first six orientations of the scheme used for the prismatic housing. In each image, rotation arrows indicate the movement required to transition from the current orientation to the following, while $x,y,z$ represent the axes of the F-frame, which in practice can be identified with the three edges formed by the intersection of the three surfaces of the 3D-printed references (Figure 2c,d). The faces of the housing are numbered to make the rotation sequence clearer.Since the IMU is rigidly attached to its housing, as shown in Figure 4, if the housing is re-oriented into N different orientations, rotation matrices ${\mathbf{R}}_i$ for $i=1,\cdots ,N$ can be computed as a combination of rotations around the axes of the F-frame. To achieve this, we begin considering that the rotation matrix ${\mathbf{R}}_i^{\prime }$, which transforms a vector expressed in the mobile frame into a vector expressed in the fixed frame when the IMU is in the i-th orientation, can be iteratively obtained from ${\mathbf{R}}_{i-1}^{\prime }$. This process depends on the axis around which the rotation is performed to transition from the $(i-1)$-th to the i-th IMU orientation:    

$${\mathbf{R}}_1^{\prime }={\mathbf{I}}_{3\times 3},{\mathbf{R}}_i^{\prime }=\left\{\begin{array}{cc}{\mathbf{R}}_x\left({\mathrm{\Delta }}_j\right){\mathbf{R}}_{i-1}^{\prime },\hfill & \mathrm{if}\mathrm{the}\mathrm{rotation}\mathrm{is}\mathrm{performed}\mathrm{around}x\hfill \\ {\mathbf{R}}_y\left({\mathrm{\Delta }}_j\right){\mathbf{R}}_{i-1}^{\prime },\hfill & \mathrm{if}\mathrm{the}\mathrm{rotation}\mathrm{is}\mathrm{performed}\mathrm{around}y\hfill \\ {\mathbf{R}}_z\left({\mathrm{\Delta }}_j\right){\mathbf{R}}_{i-1}^{\prime },\hfill & \mathrm{if}\mathrm{the}\mathrm{rotation}\mathrm{is}\mathrm{performed}\mathrm{around}z\hfill \end{array}\right.,$$

(46)

where ${\mathbf{R}}_x\left({\mathrm{\Delta }}_j\right),{\mathbf{R}}_y\left({\mathrm{\Delta }}_j\right),{\mathbf{R}}_z\left({\mathrm{\Delta }}_j\right)$ are elementary rotation matrices of angle ${\mathrm{\Delta }}_j$ around $x^{\prime },y^{\prime }$ and $z^{\prime }$ axis, respectively, while ${\mathrm{\Delta }}_j,j=i-1$ is the j-th rotation angle, required to transition from the $(i-1)$-th orientation to the i-th. Then, since Equations (7) and (26) require the acceleration and the angular velocity to be expressed in the moving frame, we obtain the needed rotation matrices as ${\mathbf{R}}_i={\mathbf{R}}_i^{\prime T}$. The computation of ${\mathbf{R}}_i$ matrices for the re-orientation sequences used in this paper is detailed in Appendix B.

The prismatic housing enables a simpler re-orientation procedure. However, if a reference non-sloped w.r.t. $\mathbf{g}$ is used (e.g., the one in Figure 2c), this re-orientation scheme is hardly optimal for accelerometer calibration. In fact, if the reference rests on a surface reasonably orthogonal to $\mathbf{g}$, one of the sensor axes would be almost parallel to $\mathbf{g}$, and the others almost orthogonal: the measures obtained on each channel would be whether near the maximum possible or near zero, meaning that they would be characterized by very small SNR. To overcome this problem, one can: (i) use a sloped reference, so that all axes measurements have greater SNR, or (ii) use a housing characterized by more possible orientations, improving also gyroscope calibration by performing more rotations. Thus, experiments on all four possible combinations of reference surface and housing were conducted. Figure 5 illustrates the experimental equipment in the initial orientation of the procedure for the Sloped 18-faced setup. A flowchart summarizing the steps needed to perform IMU calibration using our method can be found in Appendix C.

The Adafruit Precision NXP 9-DOF Breakout Board, which is a single-board IMU featuring NXP sensors FXOS8700CQ and FXAS21002, was employed for the tests. The main technical specifications of the IMU are listed in

Table 1. Main IMU technical characteristics.

Characteristic	Value
Board’s dimensions [mm]	20 × 28 × 2
Supply voltage [V]	from $1.95$ to $3.6$
Supported digital interface	$I^{2}C$ or $SPI$
Current consumption [mA]	$2.94$
Output data rate [Hz]	up to 800

, while additional details are available in [25,26,27]. The FXOS8700CQ integrates a three-axis magnetometer (which is not of interest for this work) and a three-axis digital output accelerometer, while the FXAS21002 is a three-axis digital output gyroscope. IMU’s outputs were sampled at 300 Hz to acquire raw data from the accelerometer and the gyroscope. For the gyroscope, the signal was digitally low-pass filtered with a cutoff frequency of 128 Hz before the acquisition, using the filter integrated in the sensor [26]. During the re-orientation procedure, the IMU alternates between static and non-static states: for accelerometer calibration, only raw data collected during static states were retained, while for gyroscope calibration only raw data from non-static states were kept. Non-static states were identified using a non-static detector based on the low-pass filtered gyroscope signal: to obtain this signal, the acquired one was low-pass filtered again with a cutoff frequency that, in this case, was set to $0.5$ Hz. Then, a certain time instant t was determined to be non-static if the magnitude of the filtered gyroscope signal at that instant exceeded a certain threshold, which was empirically set to 1000 Least Significant Bits (LSB). In the context of digital output sensors, an LSB represents the smallest measurable change in the physical quantity: accordingly, the physical value associated with an LSB depends on the sensor’s Full-Scale range and on the resolution $\rho =2^{\nu }$ of its Analog-to-Digital Converter (ADC), being $\nu$ the number of bits available for the ADC. For the sensor used in our tests, $1LSB$ corresponds to an angular velocity magnitude measurement of $0.01526°/\mathrm{s}$.

For the accelerometer, the mean values of the identified parameters, obtained from twelve different acquired datasets (three for each experimental setup) are reported in

Table 2. Accelerometer identification results. For each experimental setup, the left column contains mean values of identified parameters (${\theta }_l$), while the right column contains their mean percentage standard deviations (${\sigma }_{\%,{\theta }_l}$).

	Non-Sloped Prismatic		Sloped Prismatic		Non-Sloped 18-Faced		Sloped 18-Faced
	${\mathit{\theta }}_{\mathit{l}}$	${\mathbf{\sigma }}_{\%,{\mathit{\theta }}_{\mathit{l}}}$	${\mathit{\theta }}_{\mathit{l}}$	${\mathbf{\sigma }}_{\%,{\mathit{\theta }}_{\mathit{l}}}$	${\mathit{\theta }}_{\mathit{l}}$	${\mathbf{\sigma }}_{\%,{\mathit{\theta }}_{\mathit{l}}}$	${\mathit{\theta }}_{\mathit{l}}$	${\mathbf{\sigma }}_{\%,{\mathit{\theta }}_{\mathit{l}}}$
$A_{11}\left[{\displaystyle \frac{g}{LSB}}\right]$	−2.417 × ${10}^{-4}$	0.135	−2.428 × ${10}^{-4}$	0.178	−2.420 × ${10}^{-4}$	0.096	−2.422 × ${10}^{-4}$	0.257
$A_{12}\left[{\displaystyle \frac{g}{LSB}}\right]$	1.592 × ${10}^{-6}$	45.823	1.329 × ${10}^{-6}$	38.889	−4.397 × ${10}^{-6}$	3.872	−5.301 × ${10}^{-6}$	7.387
$A_{13}\left[{\displaystyle \frac{g}{LSB}}\right]$	7.193 × ${10}^{-7}$	32.321	6.985 × ${10}^{-7}$	39.377	−4.073 × ${10}^{-7}$	17.592	−1.299 × ${10}^{-6}$	30.965
$A_{22}\left[{\displaystyle \frac{g}{LSB}}\right]$	2.464 × ${10}^{-4}$	0.135	2.472 × ${10}^{-4}$	0.178	2.460 × ${10}^{-4}$	0.098	2.468 × ${10}^{-4}$	0.260
$A_{23}\left[{\displaystyle \frac{g}{LSB}}\right]$	−4.255 × ${10}^{-7}$	63.772	−4.585 × ${10}^{-7}$	60.916	−4.667 × ${10}^{-7}$	36.024	−6.040 × ${10}^{-7}$	63.096
$A_{33}\left[{\displaystyle \frac{g}{LSB}}\right]$	2.428 × ${10}^{-4}$	0.135	2.434 × ${10}^{-4}$	0.178	2.433 × ${10}^{-4}$	0.096	2.437 × ${10}^{-4}$	0.257
$b_1\left[g\right]$	2.222 × ${10}^{-3}$	35.047	2.786 × ${10}^{-3}$	36.149	1.792 × ${10}^{-3}$	40.823	1.297 × ${10}^{-3}$	101.695
$b_2\left[g\right]$	4.416 × ${10}^{-3}$	17.751	1.872 × ${10}^{-3}$	88.593	3.969 × ${10}^{-3}$	14.086	2.473 × ${10}^{-3}$	53.338
$b_3\left[g\right]$	−1.839 × ${10}^{-2}$	4.237	−1.614 × ${10}^{-2}$	5.793	−1.782 × ${10}^{-2}$	3.497	−1.474 × ${10}^{-2}$	10.020
$n_1\left[g\right]$	2.420 × ${10}^{-2}$	3.297	4.205 × ${10}^{-1}$	0.252	2.731 × ${10}^{-2}$	2.332	4.297 × ${10}^{-1}$	0.357
$n_2\left[g\right]$	−7.073 × ${10}^{-3}$	12.206	4.887 × ${10}^{-1}$	0.224	9.044 × ${10}^{-3}$	3.662	4.811 × ${10}^{-1}$	0.330
$n_3\left[g\right]$	−9.997 × ${10}^{-1}$		−7.644 × ${10}^{-1}$		−9.994 × ${10}^{-1}$		−7.641 × ${10}^{-1}$

, along with their respective mean relative standard deviations. Elements of matrix $\mathbf{A}$ are expressed in $[g/LSB]$. This measurement unit is physically meaningful for sensitivity, as it indicates how much the measured physical quantity changes in relation to the smallest possible variation in the sensor’s raw output (which, for digital sensors, is measured in LSB), while those of $\mathbf{b}$ and ${\left\{\mathbf{n}\right\}}_{M,1}$ are in $\left[g\right]$, i.e., they are normalized with respect to the modulus of the gravitational acceleration. The relative standard deviations are expressed as percentages of the parameter value and are obtained as described in Section 2.1. For the gyroscope, the mean values of the identified parameters, obtained from another twelve datasets acquired simultaneously with the accelerometer ones, are reported in

Table 3. Gyroscope identification results. For each experimental setup, the left column contains mean values of identified parameters (${\vartheta }_q$), while the right column contains their mean percentage standard deviations (${\sigma }_{\%,{\vartheta }_q}$).

	Non-Sloped Prismatic		Sloped Prismatic		Non-Sloped 18-Faced		Sloped 18-Faced
	${\mathit{\vartheta }}_{\mathit{q}}$	${\mathbf{\sigma }}_{\%,{\mathit{\vartheta }}_{\mathit{q}}}$	${\mathit{\vartheta }}_{\mathit{q}}$	${\mathbf{\sigma }}_{\%,{\mathit{\vartheta }}_{\mathit{q}}}$	${\mathit{\vartheta }}_{\mathit{q}}$	${\mathbf{\sigma }}_{\%,{\mathit{\vartheta }}_{\mathit{q}}}$	${\mathit{\vartheta }}_{\mathit{q}}$	${\mathbf{\sigma }}_{\%,{\mathit{\vartheta }}_{\mathit{q}}}$
$G_{11}\left[{\displaystyle \frac{°/s}{LSB}}\right]$	1.537 × ${10}^{-2}$	0.295	1.534 × ${10}^{-2}$	0.377	1.531 × ${10}^{-2}$	0.225	1.533 × ${10}^{-2}$	0.273
$G_{12}\left[{\displaystyle \frac{°/s}{LSB}}\right]$	−3.972 × ${10}^{-4}$	8.242	−3.560 × ${10}^{-4}$	11.720	3.950 × ${10}^{-5}$	68.006	4.276 × ${10}^{-5}$	72.336
$G_{13}\left[{\displaystyle \frac{°/s}{LSB}}\right]$	−5.796 × ${10}^{-7}$	2721.260	1.850 × ${10}^{-5}$	156.690	1.067 × ${10}^{-4}$	57.991	1.792 × ${10}^{-4}$	36.388
$G_{22}\left[{\displaystyle \frac{°/s}{LSB}}\right]$	−1.571 × ${10}^{-2}$	0.294	−1.570 × ${10}^{-2}$	0.377	−1.563 × ${10}^{-2}$	0.223	−1.568 × ${10}^{-2}$	0.271
$G_{23}\left[{\displaystyle \frac{°/s}{LSB}}\right]$	2.222 × ${10}^{-5}$	150.875	1.477 × ${10}^{-5}$	635.457	9.063 × ${10}^{-6}$	520.768	1.589 × ${10}^{-5}$	133.8609
$G_{33}\left[{\displaystyle \frac{°/s}{LSB}}\right]$	−1.554 × ${10}^{-2}$	0.265	−1.557 × ${10}^{-2}$	0.339	−1.557 × ${10}^{-2}$	0.218	−1.558 × ${10}^{-2}$	0.264
$d_1[°/s]$	−9.415 × ${10}^{-1}$	9.536	−9.416 × ${10}^{-1}$	12.310	−1.192	9.177	−8.761 × ${10}^{-1}$	14.037
$d_2[°/s]$	8.196 × ${10}^{-1}$	10.851	8.990 × ${10}^{-1}$	12.673	9.468 × ${10}^{-1}$	7.452	1.177	7.757
$d_3[°/s]$	7.847 × ${10}^{-1}$	12.267	8.037 × ${10}^{-1}$	14.296	4.282 × ${10}^{-1}$	94.752	7.564 × ${10}^{-1}$	69.528
$\varphi e_3\left[rad\right]$	−5.641 × ${10}^{-3}$	66.938	−3.837 × ${10}^{-3}$	185.509	1.811 × ${10}^{-3}$	199.231	4.254 × ${10}^{-3}$	153.948
$\varphi e_2\left[rad\right]$	4.393 × ${10}^{-3}$	40.604	3.234 × ${10}^{-3}$	67.655	−1.322 × ${10}^{-3}$	132.661	−6.115 × ${10}^{-4}$	213.256
$-\varphi e_1\left[rad\right]$	−1.655 × ${10}^{-3}$	117.084	−1.414 × ${10}^{-3}$	191.682	1.144 × ${10}^{-3}$	128.592	1.866 × ${10}^{-5}$	603.486

, along with their respective mean relative standard deviations. Elements of matrix $\mathbf{G}$ are expressed in $\left[{\displaystyle \frac{°/\mathrm{s}}{LSB}}\right]$, those of $\mathbf{d}$ are in $[°/\mathrm{s}]$, $\varphi$ is in $\left[\mathrm{rad}\right]$ and elements of $\mathbf{e}$ are dimensionless. The relative standard deviations are expressed as percentages of the parameter value and are calculated as described in Section 2.2.

Table 2 and Table 3 indicate that cross-sensitivities are negligible compared to sensitivities for both sensors, meaning that the sensors’ axes are well-aligned. Consequently, really high relative standard deviations are obtained for these parameters (a parameter is considered well-identified if its percentage standard deviation is below $5\%$); namely, they are non-essential parameters for the specific sensors tested. Bias terms also exhibit high percentage standard deviations, even if at first sight they appear to be non-negligible w.r.t. sensitivities. However, their impact on the mathematical model output is minimal due to their small values. In fact, according to Equation (21), elements of matrix $\mathbf{G}$ are to be multiplied by the raw gyroscope output $\mathbf{r}$, which is in the range $[-2^{\nu -1},+2^{\nu -1}]$: $\nu$ has a value that is between $\nu =12$ and $\nu =16$, meaning that $\mathbf{d}$ becomes negligible w.r.t.$\mathbf{G}\mathbf{r}$ if its components are the one in Table 3. Similar considerations apply to accelerometer bias terms. Accordingly, bias terms can be deemed non-essential parameters for both sensors. Furthermore, in gyroscope calibration, also bias drift is to be taken into account: however, this paper primarily focuses on identifying $\mathbf{G}$, since only non-static measurements are employed. Variations of the procedure that are able to account also for the bias drift effect and properly identify static gyroscope bias could be the focus of future works. Finally, $\varphi e_3,\varphi e_2,-\varphi e_1$ also exhibit a high percentage of standard deviations, but this is acceptable since their nominal values are very small. In fact, as we assumed in Section 2.2, $\varphi \approx 0$ and this can be confirmed by calculating $\varphi$ as in Equation (42), using parameters in Table 3: in the worst case (i.e., when $\varphi$ assumes the highest value) we obtain $\varphi =0.0073\mathrm{rad}=0.4205°$.

Parameters identification best practice suggests excluding non-essential parameters from the model and performing a second identification with the simplified model, which does not include them [23]. The results for the identification of essential parameters are presented in

Table 4. Accelerometer identification results for essential parameters. For each experimental setup, the left column contains mean values of identified parameters (${\theta }_l$), while the right column contains their mean percentage standard deviations (${\sigma }_{\%,{\theta }_l}$).

	Non-Sloped Prismatic		Sloped Prismatic		Non-Sloped 18-Faced		Sloped 18-Faced
	${\mathit{\theta }}_{\mathit{l}}$	${\mathit{\sigma }}_{\%,{\mathit{\theta }}_{\mathit{l}}}$	${\mathit{\theta }}_{\mathit{l}}$	${\mathit{\sigma }}_{\%,{\mathit{\theta }}_{\mathit{l}}}$	${\mathit{\theta }}_{\mathit{l}}$	${\mathit{\sigma }}_{\%,{\mathit{\theta }}_{\mathit{l}}}$	${\mathit{\theta }}_{\mathit{l}}$	${\mathit{\sigma }}_{\%,{\mathit{\theta }}_{\mathit{l}}}$
$A_{11}\left[{\displaystyle \frac{g}{LSB}}\right]$	−2.416 × ${10}^{-4}$	0.721	−2.428 × ${10}^{-4}$	0.504	−2.418 × ${10}^{-4}$	0.495	−2.422 × ${10}^{-4}$	0.472
$A_{22}\left[{\displaystyle \frac{g}{LSB}}\right]$	2.463 × ${10}^{-4}$	0.721	2.470 × ${10}^{-4}$	0.504	2.458 × ${10}^{-4}$	0.505	2.468 × ${10}^{-4}$	0.473
$A_{33}\left[{\displaystyle \frac{g}{LSB}}\right]$	2.425 × ${10}^{-4}$	0.721	2.433 × ${10}^{-4}$	0.504	2.431 × ${10}^{-4}$	0.495	2.437 × ${10}^{-4}$	0.473
$n_1\left[g\right]$			4.265 × ${10}^{-1}$	0.702			4.296 × ${10}^{-1}$	0.657
$n_2\left[g\right]$			4.824 × ${10}^{-1}$	0.641			4.813 × ${10}^{-1}$	0.605
$n_3\left[g\right]$	−1.000		−7.651 × ${10}^{-1}$		−1.000		−7.641 × ${10}^{-1}$

for the accelerometer and in

Table 5. Gyroscope identification results for essential parameters. For each experimental setup, the left column contains mean values of identified parameters (${\vartheta }_q$), while the right column contains their mean percentage standard deviations (${\sigma }_{\%,{\vartheta }_q}$).

	Non-Sloped Prismatic		Sloped Prismatic		Non-Sloped 18-Faced		Sloped 18-Faced
	${\mathit{\vartheta }}_{\mathit{q}}$	${\mathbf{\sigma }}_{\%,{\mathit{\vartheta }}_{\mathit{q}}}$	${\mathit{\vartheta }}_{\mathit{q}}$	${\mathbf{\sigma }}_{\%,{\mathit{\vartheta }}_{\mathit{q}}}$	${\mathit{\vartheta }}_{\mathit{q}}$	${\mathbf{\sigma }}_{\%,{\mathit{\vartheta }}_{\mathit{q}}}$	${\mathit{\vartheta }}_{\mathit{q}}$	${\mathbf{\sigma }}_{\%,{\mathit{\vartheta }}_{\mathit{q}}}$
$G_{11}\left[{\displaystyle \frac{°/s}{LSB}}\right]$	1.532 × ${10}^{-2}$	0.799	1.532 × ${10}^{-2}$	0.713	1.533 × ${10}^{-2}$	0.396	1.535 × ${10}^{-2}$	0.384
$G_{22}\left[{\displaystyle \frac{°/s}{LSB}}\right]$	−1.572 × ${10}^{-2}$	0.800	−1.573 × ${10}^{-2}$	0.714	−1.566 × ${10}^{-2}$	0.394	−1.571 × ${10}^{-2}$	0.382
$G_{33}\left[{\displaystyle \frac{°/s}{LSB}}\right]$	−1.545 × ${10}^{-2}$	0.705	−1.550 × ${10}^{-2}$	0.629	−1.554 × ${10}^{-2}$	0.383	−1.553 × ${10}^{-2}$	0.370

for the gyroscope. In general, all the essential parameters result in being well-identified for both sensors, since their percentage standard deviations fall below the $5\%$ threshold. As expected, datasets involving the 18-faced housing provide slightly better results for gyroscope calibration, with the identified parameters showing lower percentage standard deviations, while for the accelerometer, if the same housing is employed, experiments involving the sloped reference yield slightly better results in terms of percentage standard deviations of the estimated parameters.

In

Table 6. Comparison between sensors’ model parameters obtained with different methods and the related RMS error associated with the norm of the accelerometer’s model output. (The sign of the parameters has been omitted to facilitate comparison between the parameters obtained through different methods).

	MAN	TLS	NLS
$A_{11}\left[{\displaystyle \frac{g}{LSB}}\right]$	2.44 × ${10}^{-4}$	2.418 × ${10}^{-4}$	2.420 × ${10}^{-4}$
$A_{22}\left[{\displaystyle \frac{g}{LSB}}\right]$	2.44 × ${10}^{-4}$	2.458 × ${10}^{-4}$	2.468 × ${10}^{-4}$
$A_{33}\left[{\displaystyle \frac{g}{LSB}}\right]$	2.44 × ${10}^{-4}$	2.431 × ${10}^{-4}$	2.429 × ${10}^{-4}$
$G_{11}\left[{\displaystyle \frac{°/s}{LSB}}\right]$	1.5625 × ${10}^{-2}$	1.533 × ${10}^{-2}$	1.544 × ${10}^{-2}$
$G_{22}\left[{\displaystyle \frac{°/s}{LSB}}\right]$	1.5625 × ${10}^{-2}$	1.566 × ${10}^{-2}$	1.569 × ${10}^{-2}$
$G_{33}\left[{\displaystyle \frac{°/s}{LSB}}\right]$	1.5625 × ${10}^{-2}$	1.554 × ${10}^{-2}$	1.559 × ${10}^{-2}$
${\xi }_{RMS}\left[g\right]$	0.0132	0.0120	0.0119

, we present accelerometer and gyroscope essential parameters ($A_{pp}$ and $G_{pp}$, respectively) derived in three different manners: (i) from sensors datasheets [26,27] (MAN), (ii) using our method (TLS) and (iii) using the NLS-based method outlined in [12] (NLS). Although our method allows for a straightforward determination of the sign of sensors’ parameters, all the parameters identified by our method are reported as positive in Table 6 to facilitate comparison with manufacturer-provided values and those derived using the NLS-based method. The NLS-based method, relying on Equations (1) and (2), minimizes a cost function that depends only on the norm of the sensors’ model outputs, making it insensitive to the sign of the parameters. The results show that our method yields parameters comparable to those obtained using the NLS-based approach. Additionally, since a known physical reference is available for the accelerometer (i.e., $\parallel \mathbf{g}\parallel =1$), we used one of the calibration datasets to compute the RMS error between the norm of the gravity vector, estimated via the sensor’s model, and the known reference as follows:

$${\xi }_{RMS}=\sqrt{{\displaystyle \frac{1}{N}}\sum _{i=1}^N(\parallel {\tilde{\mathbf{g}}}_i{\parallel -1)}^2},$$

(47)

where ${\tilde{\mathbf{g}}}_i$ is the gravity vector estimation obtained by employing in Equation (9) the sensor’s raw outputs acquired in the i-th orientation of the procedure for the Non-Sloped 18-faced setup. ${\xi }_{RMS}$ values computed using the above-mentioned three different parameter sets are reported in Table 6, which demonstrates a relevant improvement in terms of ${\xi }_{RMS}$ when transitioning from the manufacturer-parameters-computed one to the one obtained using our method’s parameters. The NLS-based method appears to yield better results in terms of ${\xi }_{RMS}$; however, it is important to note that its performance heavily depends on the quality of the initial guess provided to the algorithm, which in this case was set equal to the parameters specified in the sensor’s datasheet.

Finally, cross-validations of the identification results were performed. We cross-validated the parameters obtained with the Non-sloped 18-faced setup (Table 4 and Table 5) using datasets acquired during experiments involving the other three setups. For each dataset, we computed $\mathbf{M}$ according to Equations (11) and (13), and $\mathbf{L}$ as in Equations (32) and (35)–(37). The corresponding model error was then calculated as ${ϵ}_a=\mathbf{M}\mathit{\theta }$ for the accelerometer and as ${ϵ}_g=\mathbf{L}\mathit{\vartheta }$ for the gyroscope. For datasets involving the non-sloped reference, the last element of $\mathit{\theta }$ was set to $-1$, while for those involving the sloped reference, the last three elements of $\mathit{\theta }$ were set to ${[n_1,n_2,n_3]}^T$, since only in the latter case all these three parameters result to be essential. Finally, we computed the cross-validation standard deviations as [22]:

$${\sigma }_a^{\prime }=\sqrt{{\displaystyle \frac{{ϵ}_a^T{ϵ}_a}{3N}}}$$

(48)

and

$${\sigma }_g^{\prime }=\sqrt{{\displaystyle \frac{{ϵ}_g^T{ϵ}_g}{3(N-1)}}}$$

(49)

for accelerometer and gyroscope datasets, respectively. Additionally, we also evaluated cross-validation standard deviations using manufacturer-provided sensor parameters. As manufacturer parameters, we used the ones in Table 6. Cross-validation results obtained using parameters identified by our method (${\sigma }_{a,TLS}^{\prime }$ and ${\sigma }_{g,TLS}^{\prime }$) are compared with the ones obtained using manufacturer-provided parameters (${\sigma }_{a,MAN}^{\prime }$ and ${\sigma }_{g,MAN}^{\prime }$) in

Table 7. Accelerometer and gyroscope cross-validation standard deviations for the three datasets setups.

	${\mathit{\sigma }}_{\mathit{a},\mathit{MAN}}^{\prime }\left[\mathit{g}\right]$	${\mathit{\sigma }}_{\mathit{a},\mathit{TLS}}^{\prime }\left[\mathit{g}\right]$	${\mathit{\sigma }}_{\mathit{g},\mathit{MAN}}^{\prime }[°]$	${\mathit{\sigma }}_{\mathit{g},\mathit{TLS}}^{\prime }[°]$
Non-sloped cubic	0.6699	0.0194	86.9227	2.2519
Sloped cubic	0.6684	0.0105	86.7385	1.5303
Non-sloped 18-faced	0.6740	0.0217	78.6770	1.5243

. These results clearly demonstrate that our method is an effective choice for improving accelerometer and gyroscope accuracy with no significant additional cost: in fact, using the parameters identified by our method, we were able to reduce cross-validation standard deviations by nearly two orders of magnitude compared to the case where manufacturer parameters were employed.

4. Conclusions

In this paper, we proposed an identification method for accelerometers and gyroscopes, namely a general identification method for six-axis IMUs. Our method can be performed in-field by non-expert users, requiring only simple 3D-printed equipment. It is based on a pre-determined sequence of rotations and the solution of two TLS problems. No prior knowledge regarding references orientation w.r.t. $\mathbf{g}$ or regarding rotation axes orientation w.r.t. gyroscope sensitive triad is required. Our identification method does not require prior knowledge of the parameters to be determined: this makes it a valuable option to generate a reliable initial parameter set to be used for NLS-based methods when manufacturer-provided parameters are unavailable (e.g., for prototypes or newly designed sensors). Experimental results showed that the proposed method is able to correctly identify sensor model parameters and that the accuracy of sensor models is severely increased w.r.t. using manufacturer-provided parameters. The proposed method is not suitable for estimating gyroscope static bias terms: this could be the focus of future works.