## Abstract

Motion reconstruction and navigation require accurate orientation estimation. Modern orientation estimation methods utilize filtering algorithms, such as the Kalman filter or Madgwick's algorithm. However, these methods do not address potential sensor saturation, which may occur within short time periods in highly dynamic applications, such as, e.g., particle tracking in snow avalanches, leading to inaccurate orientation estimates. In this paper, we present two algorithms for orientation estimation combining magnetometer and partially saturated gyrometer readings. One algorithm incorporates magnetic field vector observations and the full nonlinearity of the exponential map. The other, computationally more efficient algorithm builds on a linearization of the exponential map and is solved analytically. Both algorithms are then applied to measurement data from four different experiments, with two of them being snow avalanche experiments. Moreover, Madgwick's filtering algorithm was used to validate the proposed algorithms. The two algorithms improved the orientation estimation significantly in all experiments. Hence, the proposed algorithms can improve the performance of existing sensor fusion algorithms significantly.

## 1 Introduction

Orientation estimation of a rigid body utilizing inertial measurement units (IMUs) is crucial and a challenging task in the field of motion reconstruction and navigation and requires special computational methods. Limitations arise due to errors in IMU readings, i.e., erroneous sensor coordinate systems, sensor drift, saturation of sensor axes, and others [1]. These limitations are partly addressed with calibration [2,3] of the IMU, with sensor fusion methods such as the Kalman filter [4] and Madgwick's filter [5], and by incorporation of motion constraints [6–8]. However, the saturation of gyrometer readings, which occurs in fast-rotating applications and results in cutoff measurement data, was only addressed for the case of rotation about a single sensor axis recently [9]. Despite the solution for the case of one-dimensional rotation, a validated solution for the general, spatial rotation as it may occur in, e.g., particle tracking in snow avalanches [10,11] or rock falls [12], was not yet proposed.

The three individual sensors in a conventional 9-axis IMU, i.e., accelerometer, gyrometer, and magnetometer, are used to measure acceleration, angular velocity, and magnetic flux density, respectively. The latter sensors inherit limits regarding the maximum measurement value. Motion above these limits results in sensor saturation and hence cutoff measurement values. Magnetometer saturation can be neglected as the Earth's magnetic field strength is approximately constant and far from the saturation limit [13], whereas saturation of the accelerometer and gyrometer is possible in highly dynamic applications. Winkler et al. [14] addressed the issue of gyrometer saturation; however, there is no reference data to evaluate the performance of their proposed method and saturated gyrometer readings were still incorporated for orientation estimation.

Orientation estimation in the field of navigation or motion reconstruction is mostly conducted via sensor fusion. There is the Kalman filter and different variations of it, e.g., extended Kalman filter and unscented Kalman filter [15], where states are estimated by combining measurement data and a system model [16]. For orientation estimation, there is Wahba's problem [17], which seeks to find a rotation matrix utilizing vector observations to perform a coordinate transformation. Solutions to Wahba's problem are, e.g., QUEST [18], ESOQ [19], FOAM [20], and OLAE [21]. A popular algorithm for orientation estimation, especially from IMU readings, is Madgwick's filtering algorithm [5,22]. However, neither of these algorithms addresses gyrometer saturation, resulting in insufficient orientation estimation [23].

In this paper, we propose an algorithm that incorporates magnetometer readings to recover saturated gyrometer readings and thus the orientation of IMUs. This stand-alone algorithm can also be applied to state-of-the-art sensor fusion algorithms to improve their accuracy. To evaluate the proposed algorithm, we utilize four different data sets, one generated (artificial) IMU dataset, a dataset from a field experiment, and data from two snow avalanches. To further validate the performance of the proposed algorithm, we provide results obtained with Madgwick's filter. The remainder of this paper is structured as follows: In Sec. 2, we derive the methods for orientation estimation of gyrometer readings and—in case of saturation—magnetometer readings. Section 3 addresses the acquisition of artificial and real IMU readings. The results of the proposed algorithm, also applied to sensor fusion algorithms, are discussed in Sec. 4.

## 2 Orientation Estimation

The computation of orientation from IMU measurement data with occasionally saturated gyrometers requires the derivation of two algorithms. One, for the case of no saturation, in which orientation will be computed utilizing a Lie-group approach. Another is the case of gyrometer saturation, in which the missing information about orientation is provided by magnetometer data.

This section is structured as follows: In Sec. 2.1, we introduce the sensor data, corresponding frames, and coordinate transformations. Section 2.2 includes the algorithms to estimate orientation from angular velocity. In Sec. 2.3, we propose an algorithm for orientation estimation from magnetic flux densities and partly saturated angular velocities. Section 2.4 relates the presented approach to Wahba's problem.

### 2.1 Data Representation.

*t*and $ti+1$ denote the

_{i}*i*-th and (

*i*+ 1)-th time-step, respectively, where $i\u2208{0,1,2,\u2026,n\u22121}$. The stepsize $\Delta \u2009t$ equals 2.5 ms with an accuracy of $3\u2009\mu s$ for the IMU data in this paper. In addition, IMU data are represented in different frames (denoted by a left superscript), as the IMU is moving over time. There is the initial frame (I) corresponding to time-step

*t*

_{0}and the sensor fixed frame (S). For a coordinate transformation from the sensor frame (S) to the initial frame (I) for the

*i*-th time-step, we introduce the transformation

**R**$\u2208$ SO(3)

and $ri$ being a placeholder representing a spatial vector.

### 2.2 Orientation From Angular Velocity.

### 2.3 Orientation From Magnetic Flux Density.

^{2}

*z*-component of the

*i*-th time-step is saturated, we denote the incremental rotation vector

To estimate the unknown components of the incremental rotation vector $S\Omega i\u2032$ and hence angular velocity $S\omega i$, we derive two methods; one incorporating the nonlinearity of the exponential map and another with a linearized exponential map, defined in Secs. 2.3.1 and 2.3.2, respectively.

#### 2.3.1 Nonlinear Estimation.

*z*-axis, fsolve is initialized with a guess for the saturated incremental rotation

where $S\omega z,i\u22121$ is the *z*-component of the angular velocity from the previous time-step.

#### 2.3.2 Linear Estimation.

A computationally more efficient variant of the method presented in Sec. 2.3.1 can be derived by linearization of the exponential map (Eq. (7)) and is shown in this section.

*z*-axis, we obtain two equations from Eq. (17) to estimate the

*z*-component of $S\Omega \u2032$

Experiments on generated data show that the accuracy of the estimation depends on the magnitude of the denominator from Eqs. (18) and (19), respectively $Shy,i+1$ and $Shx,i+1$, which is shown in Fig. 1. Higher accuracy is achieved by using the equation with larger denominator value, which is therefore implemented in the present paper. Note that there is a limitation of the presented estimation algorithm if the magnetic flux density vector $Sh$ is parallel to the saturated angular velocity component, as seen in Eqs. (18) and (19) for the *z*-component.

Finally, a scheme of the algorithm for orientation estimation is displayed in Fig. 2, showing the different equations used, whether or not, the gyrometer is saturated.

### 2.4 Relation to Wahba's Problem.

**R**$\u2208$ SO(3) to transform between two coordinate systems through minimization of a loss function in the form of [31]

where $bk,i$ and $bk,i+1$ are two vectors corresponding to different time steps *i* and (*i *+* *1), respectively, and *w _{k}* is a term to weigh the

*k*individual vector observations. Wahba's problem is of particular interest for the present work as the presented algorithms in Sec. 2.3 also seek for a rotation matrix $exp\u2009(S\Omega i+1)$ by observation of the vectors $Shi$ and $Shi+1$, see Eq. (12). However, as there have to be at least two vector observations, i.e., $N\u22652$, the present solutions addressing Wahba's problem, e.g., QUEST, ESOQ, FOAM, and OLAE are not sufficient in the case of saturated gyrometer data. This is due to the fact that the measured accelerations are comprised of gravitation and translational acceleration and separation of the latter two components depends on $Ri$ and thus $\omega i$ [6], see Eqs. (5) and (6). Hence, there is only one vector observation left, which is the magnetic flux density measurement.

## 3 IMU Data Acquisition

To evaluate the performance of the proposed algorithm, four IMU datasets were acquired. One with IMU data computed within a rigid body simulation, simplifying real IMU data regarding stochastic and deterministic sensor errors. In addition, an experiment utilizing a rotating platform under known conditions was performed and in situ measurements from two snow avalanches were acquired. The latter three experiments were measured with the AvaNode [32], which uses an IMU of type MPU9250 with saturation limits of ±16 × gravity, ±2000 deg/s, and ±4.8 mT, for the accelerometer, gyrometer, and magnetometer, respectively.

Note that for the rotating platform experiment, we performed a separate experiment for magnetometer calibration, whereas for the snow avalanche experiments, we used the avalanche data itself to perform the magnetometer calibration. Gyrometers were partly calibrated incorporating the bias term from Neurauter et al. [6] but neglecting the calibration matrix. Accelerometers were not calibrated as translational accelerations are only used for comparison purposes with Madgwick's algorithm.

To validate the angular velocity estimates of the proposed algorithms, we artificially saturated the gyrometer readings and used the actual, unsaturated values as a reference. Saturation results in data that are cut off above the saturation limit such that every sample in saturation takes the value of the saturation limit.

### 3.1 Simulated Measurement Data.

where $Ri$ are the transformations derived from Eq. (5). Due to simulated data, constant Euclidean norm of the magnetic flux density vector (Eq. (9)) is satisfied, see Fig. 4. Note that we removed the first 0.14 s of the measurement such that the measurement does not start in saturation which would cause larger errors for the first estimate as there is no sufficient angular velocity from the previous time-step, see Eq. (15).

### 3.2 Rotating Platform Data.

To test the algorithm on real IMU data, a field experiment was conducted. Here, we used a rotating platform to rotate the AvaNode about a single axis without translational movement, see Fig. 5. As magnetic field measurements are sensitive to metals, the experimental setup was built from wood and 3D-printed plastic parts. In addition, only the IMU was placed on the rotating table and the rest of the AvaNode (electronics, batteries, etc.) was placed under the rotating table such that the latter two were approx. 0.5 m apart. Prior to the measurement, the IMU was aligned such that the *y*-axis pointed toward magnetic north while the *z*-axis pointed upwards vertically. These idealizations and calibration of the magnetometer at the experimental site yield magnetic flux densities with an approximately constant Euclidean norm as shown in Fig. 6. In Fig. 7, the measured angular velocities from the rotating platform experiment are shown with artificial saturation limits at ±5 rad/s.

### 3.3 Snow Avalanche Data.

This section covers two snow avalanches and corresponding datasets, which were acquired with the AvaNode. In the first snow avalanche experiment (Innsbruck), no gyrometer saturation occurred and thus we saturated artificially. In the second snow avalanche experiment (Wattener Lizum), actual gyrometer saturation occurred.

#### 3.3.1 Snow Avalanche 1: Innsbruck.

This experiment took place on Mar. 15, 2021, at the Nordkette in Innsbruck, Austria, at a location denoted as Seilbahnrinne (47 deg18′44″ N, 11 deg22′60″ E, 2269 m a.s.l.) shown in Fig. 8. Here, the AvaNode was dropped from the avalanche control tower into the release area before the snow avalanche was artificially triggered by avalanche control measures, descended, and stopped at the deposition area shown in Fig. 8. This snow avalanche can be split into an acceleration phase (2 s–6 s), a steady-state flow phase (6 s–35 s), and the deceleration phase (34 s–45 s), see Fig. 9. As a first step toward angular velocity reconstruction in snow avalanches, we consider the acceleration phase. Considered angular velocities of the acceleration phase are displayed in Fig. 10. Corresponding magnetic flux densities are shown in Fig. 11.

#### 3.3.2 Snow Avalanche 2: Wattener Lizum.

This experiment took place on Mar. 10, 2021, at the Wattener Lizum, Austria (47 deg10′24″ N, 11 deg35′59″ E, 2126 m a.s.l). Here, the AvaNode was manually placed in the potential release area, where the avalanche was artificially triggered with large additional load. Thereafter, a snow avalanche was triggered by blasting. The traveled distance of the AvaNode was approximately 100 m. Measured magnetic flux densities and angular velocities are shown in Figs. 12 and 13, respectively. In contrast to the previously presented measurements, actual gyrometer saturation occurred in this experiment, see Fig. 13. Thus, angular velocities above 2000 deg/s (34.9 rad/s) were not measured and hence, orientation cannot be determined. Although the magnetometer is calibrated, we see a significant deviation in the norm when the AvaNode is rotating fast, see Fig. 12.

## 4 Results and Discussion

In this section, we introduce the utilized implementation of Madgwick's filter in Sec. 4.1 and show the results obtained with the proposed algorithms applied to simulated IMU data in Sec. 4.2. Moreover, the results for real IMU measurement data from a rotating table experiment and two snow avalanches are shown and discussed in Secs. 4.3 and 4.4, respectively.

The proposed algorithms from Sec. 2 are applied to four acquired datasets described in Sec. 3. Evaluation of the estimation algorithms is conducted in the angular velocity domain as well as in the orientation domain. The estimated angular velocities from Eqs. (14) and (17), referred to as nonlinear and linear, respectively, are compared to the reference data. In the orientation domain, Madgwick's filtering algorithm was applied to saturated and estimated data.

### 4.1 Madgwick's Filtering Algorithm.

To evaluate the presented method for the estimation of saturated angular velocity components, we use Madgwick's filtering algorithm, which estimates the orientation of an object by combining data from an accelerometer, a gyrometer, and a magnetometer. The algorithm derives an orientation estimate by minimization (using gradient descent) of a cost function based on the difference between the estimated orientation and the measurement data from the previous time-step [5,22].

^{3}, which employs Euler parameters (unit quaternions) $q\u2208S3$ [35]

**q**. Deviating from the standard settings, we adjusted the filter gain, sampling frequency, and initial orientation according to Table 1. To underline the significance of the presented algorithms, five different methods for orientation estimation are compared, which are defined in Table 2. For evaluation of the individual methods, we compare the angular error in the orientation

*α*, which is defined as

where $Rr$ are reference rotation matrices derived from Madgwick's filter with reference angular velocities as input and $Re$ are the estimated rotation matrices, see Eq. (5). The matrix logarithm in Eq. (28) is a Euler parameter implementation in Exudyn (version 1.6.119.dev1) [36] of the matrix logarithm defined in [37].

Notation | Variable | Setting |
---|---|---|

Filter gain | — | Gain = 0.01 |

Sampling frequency | 1/$\Delta \u2009t$ | Freq = 400 |

Initial orientation | $q0$ | q0 = $[1,0,0,0]$ |

Notation | Variable | Setting |
---|---|---|

Filter gain | — | Gain = 0.01 |

Sampling frequency | 1/$\Delta \u2009t$ | Freq = 400 |

Initial orientation | $q0$ | q0 = $[1,0,0,0]$ |

A: $\omega s$ → Eq. (14) → $\Omega e$ → Eq. (5) → $Re$ |

B: $\omega s$ → Eq. (14) → $\Omega e$ → Eq. (20) → $\omega e$ → MF → $qe$ → $Re$ |

C: $\omega s$ → Eq. (17) → $\Omega e$ → Eq. (5) → $Re$ |

D: $\omega s$ → Eq. (17) → $\Omega e$ → Eq. (20) → $\omega e$ → MF → $qe$ → $Re$ |

E: $\omega s$ → MF → $qe$ → $Re$ |

A: $\omega s$ → Eq. (14) → $\Omega e$ → Eq. (5) → $Re$ |

B: $\omega s$ → Eq. (14) → $\Omega e$ → Eq. (20) → $\omega e$ → MF → $qe$ → $Re$ |

C: $\omega s$ → Eq. (17) → $\Omega e$ → Eq. (5) → $Re$ |

D: $\omega s$ → Eq. (17) → $\Omega e$ → Eq. (20) → $\omega e$ → MF → $qe$ → $Re$ |

E: $\omega s$ → MF → $qe$ → $Re$ |

Indizes s and e denote saturated and estimated, respectively. Method A denotes the nonlinear estimate. Method B denotes the linear estimate. Method C denotes MF used in combination with the nonlinear estimate. Method D denotes MF used in combination with the linear estimate. Method E denotes the original MF used in combination with saturated angular velocities.

### 4.2 Simulated Measurement Evaluation.

This section covers the estimation of angular velocity and orientation in the case of one angular velocity component being saturated, hence considering the saturation limit 1 from Fig. 3 at ±39 rad/s. Applying the methods from Secs. 2.3.1 and 2.3.2 to saturated angular velocities from Fig. 3 respectively yields the nonlinear and linear estimates shown in Fig. 14 with the corresponding errors shown in Fig. 15. The latter errors are below $10\u221211$ rad/s for the nonlinear estimate and 2 rad/s (approximately 5% with respect to the reference) for the linear estimate. These errors depend on the stepsize $\Delta \u2009t$, see Eq. (20). Additionally, errors in the linear estimate are due to the linearization of the rotation, see Eq. (16). The performances of the individual methods regarding orientation estimation are displayed in Fig. 16 with corresponding mean and maximum errors displayed in Table 3, showing that method B yields the most accurate estimate with errors below $10\u221214$. In contrast, method E provides the least accurate estimate. Figure 16 further shows that methods A and C improve the orientation estimate without further processing, reducing the mean error in orientation estimation by 94.5% and 92.4%, respectively, compared to the state-of-the-art method E. However, the application of Madgwick's filter on methods A and C further improves the orientation estimate, see methods B and D.

Method | Mean error (rad) | Max error (rad) |
---|---|---|

A | 4.205 03 × 10^{−2} | 2.037 95 × 10^{−1} |

B | 1.188 15 × 10^{−15} | 6.745 36 × 10^{−15} |

C | 5.811 95 × 10^{−2} | 2.021 03 × 10^{−1} |

D | 1.961 52 × 10^{−2} | 7.528 82 × 10^{−2} |

E | 7.699 88 × 10^{−1} | 3.058 41 |

Method | Mean error (rad) | Max error (rad) |
---|---|---|

A | 4.205 03 × 10^{−2} | 2.037 95 × 10^{−1} |

B | 1.188 15 × 10^{−15} | 6.745 36 × 10^{−15} |

C | 5.811 95 × 10^{−2} | 2.021 03 × 10^{−1} |

D | 1.961 52 × 10^{−2} | 7.528 82 × 10^{−2} |

E | 7.699 88 × 10^{−1} | 3.058 41 |

#### 4.2.1 Saturation of Two Components at a Time.

In this section, we present the results of the estimation in the case of two angular velocity components at a time being saturated, hence using the saturation limit 2 from Fig. 3 at ±30 rad/s. Note that for readability purposes, we do not show the results of the *z*-component estimation in this section as it was already addressed in Sec. 4.2.

Applying the methods from Secs. 2.3.1 and 2.3.2 to saturated angular velocities from Fig. 3 respectively yields the linear estimates shown in Fig. 17 and the errors in the nonlinear estimates shown in Fig. 18. The errors for the nonlinear estimate of components *x* and *y* are below $1\xd710\u221211$, see Fig. 18. As already shown in Fig. 1, the errors of the linear estimate (method C) are indirectly proportional to the magnitude of the magnetic field component. In the case of two saturated axes at a time, this cannot be circumvented by choosing the equation with the larger denominator, as there are two equations with one unknown respectively to solve the two unknown components of the incremental rotation vector. Hence, the linear estimate shown in Fig. 17 diverges when the magnetic flux density equals zero. However, for orientation estimation, method C still yields better results than the state-of-the-art method E, reducing the mean error by 65.6% see Table 4. Method D has no increased accuracy in the orientation estimation compared to method C in the case of two saturated axes at a time. The proposed method B is, as expected, the most accurate with errors below $10\u221212$.

Method | Mean error (rad) | Max error (rad) |
---|---|---|

A | 4.205 03 × 10^{−2} | 2.037 95 × 10^{−1} |

B | 3.717 18 × 10^{−14} | 2.906 79 × 10^{−13} |

C | 2.906 04 × 10^{−1} | 9.694 08 × 10^{−1} |

D | 2.882 28 × 10^{−1} | 9.898 53 × 10^{−1} |

E | 8.456 75 × 10^{−1} | 3.019 77 |

Method | Mean error (rad) | Max error (rad) |
---|---|---|

A | 4.205 03 × 10^{−2} | 2.037 95 × 10^{−1} |

B | 3.717 18 × 10^{−14} | 2.906 79 × 10^{−13} |

C | 2.906 04 × 10^{−1} | 9.694 08 × 10^{−1} |

D | 2.882 28 × 10^{−1} | 9.898 53 × 10^{−1} |

E | 8.456 75 × 10^{−1} | 3.019 77 |

### 4.3 Rotating Platform Evaluation.

Applying the methods from Secs. 2.3.1 and 2.3.2 to saturated angular velocities from Fig. 7, respectively, yields the nonlinear and linear estimates shown in Fig. 20 with the corresponding errors shown in Fig. 21. The errors in Fig. 21 are large compared to the errors in simulated data shown in Fig. 15, which indicates that the origin of the error can be related to the measurement data. However, the mean value of the errors from Fig. 21 equals $\u22123.4\xd710\u22123$ rad/s and $4.4\xd710\u22123$ rad/s for the nonlinear and linear angular velocity estimate, respectively. Hence, the derived terminal angles are comparably accurate with an error of $\u22121.6\xd710\u22122$ rad and $2.1\xd710\u22122$ rad for the nonlinear and linear estimate, respectively, see Fig. 22. The performances of the individual methods regarding orientation estimation are displayed in Fig. 23 with corresponding mean and maximum errors shown in Table 5. Method B yields the most accurate orientation estimation followed by method D with errors in the same range. The nonlinear and linear estimates without Madgwick's filter, hence methods A and C, respectively, also yield errors below $10\u22123$. The proposed methods A and C reduce the mean and maximum errors for orientation estimation compared to the state of the art method E by approximately 84% and 92%, respectively, see Table 5. The proposed methods B and D reduce the mean and maximum errors for orientation estimation compared to method E by approximately 95% and 96%, respectively.

Method | Mean error (rad) | Max error (rad) |
---|---|---|

A | 7.592 44 × 10^{−4} | 3.280 45 × 10^{−3} |

B | 2.239 75 × 10^{−4} | 1.855 92 × 10^{−3} |

C | 7.688 87 × 10^{−4} | 3.459 96 × 10^{−3} |

D | 2.408 54 × 10^{−4} | 1.437 39 × 10^{−3} |

E | 4.742 50 × 10^{−3} | 4.454 77 × 10^{−2} |

Method | Mean error (rad) | Max error (rad) |
---|---|---|

A | 7.592 44 × 10^{−4} | 3.280 45 × 10^{−3} |

B | 2.239 75 × 10^{−4} | 1.855 92 × 10^{−3} |

C | 7.688 87 × 10^{−4} | 3.459 96 × 10^{−3} |

D | 2.408 54 × 10^{−4} | 1.437 39 × 10^{−3} |

E | 4.742 50 × 10^{−3} | 4.454 77 × 10^{−2} |

### 4.4 Snow Avalanche Evaluation.

In this section, we apply the proposed methods to data from two snow avalanches. One snow avalanche (Innsbruck) where gyrometers were artificially saturated and another snow avalanche (Wattener Lizum) where actual gyrometer saturation occurred. Hence, for the latter snow avalanche, no reference data is available.

#### 4.4.1 Snow Avalanche 1: Innsbruck.

Applying the methods from Secs. 2.3.1 and 2.3.2 to saturated angular velocities from Fig. 10 respectively yields the nonlinear and linear estimates shown in Fig. 24 with the corresponding errors shown in Fig. 25. The maximum errors are 1.3 rad/s and 1.1 rad/s for the linear and nonlinear estimates, respectively.

The performance of the individual methods regarding orientation estimation is shown in Fig. 26 with corresponding mean and maximum errors displayed in Table 6. All methods perform significantly worse for the present avalanche data compared to the rotating table experiment, compare Figs. 26 and 23. We assume that this is due to the idealizations discussed in Sec. 3.2, through which there was no influence of the electronics from the AvaNode regarding magnetic field measurements during the rotating table experiment. However, methods A–D still reduce the mean error by approximately 90% and the maximum error by approximately 86% compared to the state-of-the-art method E, see Table 6.

#### 4.4.2 Snow Avalanche 2: Wattener Lizum.

Applying the methods from Secs. 2.3.1 and 2.3.2 to saturated angular velocities from Fig. 13 respectively yields the nonlinear and linear estimates shown in Fig. 27. As actual saturation occurred in this experiment we cannot define the error in the estimates or validate the performance of the orientation estimation, as there is no reference. However, we assume that the magnitudes of the angular velocity estimates are in a realistic range. Despite the realistic magnitudes, the sign of the estimates seems to be partially wrong, for instance, see Fig. 27 at approximately 4.5 s. This might be due to the fluctuation of the magnetic field norm, see Fig. 12.

## 5 Conclusion

Saturation of gyrometer readings and thus cutoff angular velocity data is a common issue when utilizing low-cost IMUs for orientation estimation of fast-rotating objects, e.g., particles in snow avalanches. Therefore, we proposed an algorithm to estimate the saturated angular velocities by incorporation of magnetic flux density readings and consideration of unsaturated angular velocity components. This algorithm takes advantage of the constant direction and magnitude of the magnetic field vector to derive the missing information regarding rotation by minimization of a nonlinear objective function. In addition, we simplified the algorithm deriving a linearized system of equations, which can be computed analytically. Evaluation of the algorithms was conducted utilizing four data sets from different experiments, i.e., a simulated measurement, an experiment using a rotating platform, and two snow avalanches. Moreover, Madgwick's algorithm, a popular sensor fusion algorithm, was combined with the presented algorithms. For generated data, the nonlinear algorithm estimated the saturated angular velocity components with a remaining error of $1\xd710\u221212$ rad/s and thus validating the performance of the proposed algorithms. The linear algorithm compared to the nonlinear algorithm performed significantly worse for simulated data but only slightly worse for measurement data. Hence, the main error is due to the IMU, specifically the magnetometer. Nevertheless, results derived from Madgwick's algorithm used on estimated angular velocities point out, that both, the proposed linear and nonlinear algorithms provide more accurate orientation estimation than Madgwick's algorithm used on saturated data. Moreover, the mean and maximum errors in orientation estimation of a particle in a snow avalanche were reduced by 90% and 86%, respectively, compared to the state of the art. Thus, the proposed algorithms provide reasonable angular velocity estimates and may be implemented to state of the art sensor fusion algorithms to increase the accuracy in case of gyrometer saturation. In addition, the linear algorithm could be used for real-time estimation applications.

## Funding Data

This work was conducted as part of the international cooperation project “AvaRange - Particle Tracking in Snow Avalanches” supported by the German Research Foundation (DFG, project No. 421446512) and the Austrian Science Fund (FWF, project No. I 4274-N29).

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## Notes

Equation (10) holds for arbitrary vectors, but for the presented approach we require the vector to have a constant direction to derive a solution for the rotation.

https://github.com/Mayitzin/ahrs (accessed on May 9, 2023).

## References

**41**(2), pp. 261–280.https://www.researchgate.net/publication/4718386_Attitude_Determination_from_Vector_Observations_A_Fast_Optimal_Matrix_Algorithm