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 [68]. 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.

In this paper, we utilize spatial magnetic flux densities h3 and spatial angular velocities ω3 measured by a magnetometer and gyrometer of an IMU, respectively. The sampling frequency of the IMU determines the stepsize
Δti=(ti+1ti)
(1)
where ti and ti+1 denote the i-th and (i + 1)-th time-step, respectively, where i{0,1,2,,n1}. The stepsize Δt equals 2.5 ms with an accuracy of 3μ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 t0 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
Iri=RiSri
(2)
with the rotation matrix R SO(3)
SO(3):={R3×3|RRT=I,det(R)=+1}
(3)

and ri being a placeholder representing a spatial vector.

2.2 Orientation From Angular Velocity.

To determine the orientation utilizing angular velocities measured with a gyrometer, we solve the kinematic reconstruction equation
R˙=Rω˜
(4)
where R˙ is the time derivative of R and ω˜so(3) describes the skew symmetric matrix of angular velocities ω, such that ω×y=ω˜y for ω,y3 [24]. For a solution of Eq. (4), we utilize an approach in the form of [25,26]
Ri+1=Riexp(SΩ˜i)
(5)
with the incremental rotation vector
SΩi=ΔtSωi
(6)
represented in the sensor frame (S) and the matrix exponential exp(Ω˜) using the Euler–Rodrigues formula [27]
exp(Ω˜)=I+sinc(||Ω||)Ω˜+12sinc2(||Ω||2)Ω˜2
(7)
with the cardinal sine function [28]
sinc(||Ω||)={1if||Ω||=0sin||Ω||||Ω||else
(8)

2.3 Orientation From Magnetic Flux Density.

As the Earth's magnetic field vector h3 with respect to the initial sensor frame (I) is constant for short time measurements
||Ih||2=const.
(9)
we obtain2
Ihi=RiShi
(10)
using Eq. (2) and hence [14]
Ri+1Shi+1=RiShi
(11)
Inserting Eq. (5) on the left-hand side of Eq. (11) and eliminating the rotation matrices Ri yields
exp(SΩ˜i)Shi+1=Shi
(12)
Equation (12) represents a nonlinear system of equations, where—in case of saturation—the individual angular velocities comprised in the incremental rotation vector are partially or fully unknown, which is indicated by an apostrophe, SΩ. For instance, if the z-component of the i-th time-step is saturated, we denote the incremental rotation vector
SΩi=[SΩx,i,SΩy,i,SΩz,i]T
(13)

To estimate the unknown components of the incremental rotation vector SΩi and hence angular velocity Sω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.

To derive a solution for the unknown components of the incremental rotation vector, we define the residual R
R(SΩi)=||exp(SΩ˜i)Shi+1Shi||2
(14)
where one or two components of SΩi are unknown due to saturation, and minimize it with respect to the saturated components of SΩi utilizing fsolve, an implementation of Newton's method, from scipy (version 1.7.3), a package containing fundamental algorithms for scientific computing in Python [29]. Standard settings are used. Following the example of the saturated z-axis, fsolve is initialized with a guess for the saturated incremental rotation
SΩz,i=ΔtSωz,i1
(15)

where Sωz,i1 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.

The series expansion of the exponential map reads [30]
exp(Ω˜)=I+Ω˜1!+Ω˜22!+
(16)
where we consider up to the linear term. By substitution of Eq. (16) into Eq. (14) and after elementary transformations, we obtain a system of linear equations for the unknown SΩ
Sh˜i+1SΩ=ShiShi+1
(17)
It is important to note that the coefficient matrix Sh˜i+1 of the linear system Eq. (17) does not have full rank, i.e., rank(Sh˜i+1)=2. Thus, Eq. (17) has no unique solution in case all (three) components of Sωi+1 are saturated. However, at least two components of SΩ can be reconstructed. Note, while there is one equation for each unknown component of SΩ in case two components of the angular velocity vector Sωi+1 are saturated, Eq. (17) also gives two equations in case only one component of Sωi+1 is saturated. For example, in case of a saturated z-axis, we obtain two equations from Eq. (17) to estimate the z-component of SΩ
SΩz,i=Shx,i+Shx,i+1+SΩy,iShz,i+1Shy,i+1
(18)
SΩz,i=Shy,iShy,i+1+SΩx,iShz,i+1Shx,i+1
(19)
Equations (5) and (6) only give a first-order approximation to the solution of the kinematic reconstruction equation Eq. (4), cf. [28]. Thus, the angular velocity vector Sωi reconstructed from the incremental rotation vector estimate SΩ is also a first-order approximation whose error depends linearly on the time-step size Δt, i.e.,
Sωi=SΩiΔt+O(Δt),(Δt0)
(20)

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.

Fig. 1
Relation between the magnitude of magnetic flux densities and the linear estimation of angular velocities for the case of saturation above ±39 rad/s
Fig. 1
Relation between the magnitude of magnetic flux densities and the linear estimation of angular velocities for the case of saturation above ±39 rad/s
Close modal

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.

Fig. 2
Overview of the estimation algorithm for the i-th time-step which derives orientation from either angular velocities or magnetic flux densities and partly saturated angular velocities
Fig. 2
Overview of the estimation algorithm for the i-th time-step which derives orientation from either angular velocities or magnetic flux densities and partly saturated angular velocities
Close modal

2.4 Relation to Wahba's Problem.

In the field of orientation estimation, there is the well-known problem first defined by Wahba [17], which seeks to find a rotation matrix R SO(3) to transform between two coordinate systems through minimization of a loss function in the form of [31]
L(R)=k=1Nwk||bk,iRbk,i+1||2,N2
(21)

where bk,i and bk,i+1 are two vectors corresponding to different time steps i and (i +1), respectively, and wk 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(SΩ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., N2, 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 ω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.

Simulated IMU data are very accurate compared to real IMU measurements, which suffer from stochastic and deterministic errors. The latter can be minimized with individual calibration methods. In this paper, magnetometers were calibrated according to Renaudin et al. [3] and scaled such that the Euclidean norm of the initial magnetic flux density
||h0||2=1
(22)

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.

To obtain simulated IMU data, we simulated the free rotation of a rigid body according to [33,34]. The inertia tensor with respect to the center of mass is given by
J=diag(5.2988,1.775,4.3568)
(23)
The body fixed initial angular velocity is modified (differing from Terze et al. and Holzinger et al. [33,34]) to fit typical IMU data magnitudes yielding
Sω0=[1,0,42]T
(24)
which results in an unstable rotation due to the special choice of J, see Fig. 3. The corresponding magnetometer data, shown in Fig. 4, is generated by rotating the initial magnetic flux density vector chosen as
Sh0=[1,1,1]T
(25)
via
Shi=RiTIh0
(26)

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).

Fig. 3
Generated angular velocities and saturation limits at ±30 rad/s and ±39 rad/s for saturation of two or one angular velocity components at a time, respectively
Fig. 3
Generated angular velocities and saturation limits at ±30 rad/s and ±39 rad/s for saturation of two or one angular velocity components at a time, respectively
Close modal
Fig. 4
Generated magnetic flux densities and Euclidean norm of the latter
Fig. 4
Generated magnetic flux densities and Euclidean norm of the latter
Close modal

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.

Fig. 5
Sketch of the rotating platform experiment with the IMU placed on the rotating platform, such that the z-axis is aligned with the rotation axis
Fig. 5
Sketch of the rotating platform experiment with the IMU placed on the rotating platform, such that the z-axis is aligned with the rotation axis
Close modal
Fig. 6
Calibrated magnetic flux density as well as the Euclidean norm of the magnetic flux density from the field experiment with rotation about the z-axis
Fig. 6
Calibrated magnetic flux density as well as the Euclidean norm of the magnetic flux density from the field experiment with rotation about the z-axis
Close modal
Fig. 7
Measured angular velocities from the field experiment with rotation about the z-axis and saturation limits at ±5 rad/s
Fig. 7
Measured angular velocities from the field experiment with rotation about the z-axis and saturation limits at ±5 rad/s
Close modal

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.

Fig. 8
Snow avalanche test area at the Nordkette in Innsbruck, Austria, called Seilbahnrinne with markers indicating important areas
Fig. 8
Snow avalanche test area at the Nordkette in Innsbruck, Austria, called Seilbahnrinne with markers indicating important areas
Close modal
Fig. 9
Measured angular velocities from the snow avalanche at the Nordkette showing the acceleration phase, steady-state flow phase, and deceleration phase
Fig. 9
Measured angular velocities from the snow avalanche at the Nordkette showing the acceleration phase, steady-state flow phase, and deceleration phase
Close modal
Fig. 10
Measured angular velocities from the acceleration phase of a snow avalanche and saturation limits at ±5 rad/s
Fig. 10
Measured angular velocities from the acceleration phase of a snow avalanche and saturation limits at ±5 rad/s
Close modal
Fig. 11
Calibrated magnetic flux density as well as the Euclidean norm of the magnetic flux density from the acceleration phase of a snow avalanche
Fig. 11
Calibrated magnetic flux density as well as the Euclidean norm of the magnetic flux density from the acceleration phase of a snow avalanche
Close modal

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.

Fig. 12
Calibrated magnetic flux density as well as the Euclidean norm of the magnetic flux density from the snow avalanche Wattens
Fig. 12
Calibrated magnetic flux density as well as the Euclidean norm of the magnetic flux density from the snow avalanche Wattens
Close modal
Fig. 13
Measured angular velocities from the snow avalanche in Wattens with actual gyrometer saturation
Fig. 13
Measured angular velocities from the snow avalanche in Wattens with actual gyrometer saturation
Close modal

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].

In this work, we utilize a Python implementation of Madgwick's algorithm named AHRS3, which employs Euler parameters (unit quaternions) qS3 [35]
S3:={q=[q0,q1,q2,q3]T4:||q||2=1}
(27)
to describe the orientation. In Eq. (27), q0 represents the scalar and [q1,q2,q3]T the vector part of Euler parameters 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
α=||β||2,withβ˜=log(RrReT)
(28)

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].

Table 1

Settings deviating from the standard settings for the AHRS implementation of Madgwick's filtering algorithm, where variable describes the notation used in the present paper and setting denotes the actual parameter and corresponding values set in the Python code

NotationVariableSetting
Filter gainGain = 0.01
Sampling frequency1/ΔtFreq = 400
Initial orientationq0q0 = [1,0,0,0]
NotationVariableSetting
Filter gainGain = 0.01
Sampling frequency1/ΔtFreq = 400
Initial orientationq0q0 = [1,0,0,0]
Table 2

Different methods to derive orientation estimates differing in the utilized algorithm to compute the unknown angular velocity components and whether or not, Madgwick's filter (MF) was applied

A: ωs → Eq. (14)Ωe → Eq. (5)Re
B: ωs → Eq. (14)Ωe → Eq. (20)ωe → MF → qeRe
C: ωs → Eq. (17)Ωe → Eq. (5)Re
D: ωs → Eq. (17)Ωe → Eq. (20)ωe → MF → qeRe
E: ωs → MF → qeRe
A: ωs → Eq. (14)Ωe → Eq. (5)Re
B: ωs → Eq. (14)Ωe → Eq. (20)ωe → MF → qeRe
C: ωs → Eq. (17)Ωe → Eq. (5)Re
D: ωs → Eq. (17)Ωe → Eq. (20)ωe → MF → qeRe
E: ωs → MF → qeRe

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 1011 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 Δt, 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 1014. 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.

Fig. 14
Reference data, as well as the linear and nonlinear estimate of the z-component of the angular velocity from the generated dataset
Fig. 14
Reference data, as well as the linear and nonlinear estimate of the z-component of the angular velocity from the generated dataset
Close modal
Fig. 15
Absolute error of the linear and nonlinear estimate with respect to the reference angular velocities for the generated dataset
Fig. 15
Absolute error of the linear and nonlinear estimate with respect to the reference angular velocities for the generated dataset
Close modal
Fig. 16
Comparison of different methods for orientation estimation for the generated dataset, methods A and C being the nonlinear and linear estimates, respectively. The individual methods are defined in Table 2.
Fig. 16
Comparison of different methods for orientation estimation for the generated dataset, methods A and C being the nonlinear and linear estimates, respectively. The individual methods are defined in Table 2.
Close modal
Table 3

Mean and maximum errors of the orientation estimates shown in Fig. 16 

MethodMean error (rad)Max error (rad)
A4.205 03 × 10−22.037 95 × 10−1
B1.188 15 × 10−156.745 36 × 10−15
C5.811 95 × 10−22.021 03 × 10−1
D1.961 52 × 10−27.528 82 × 10−2
E7.699 88 × 10−13.058 41
MethodMean error (rad)Max error (rad)
A4.205 03 × 10−22.037 95 × 10−1
B1.188 15 × 10−156.745 36 × 10−15
C5.811 95 × 10−22.021 03 × 10−1
D1.961 52 × 10−27.528 82 × 10−2
E7.699 88 × 10−13.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×1011, 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 1012.

Fig. 17
Reference data, as well as the linear estimate of the x- and y-components of the angular velocity from the generated dataset with saturation limit 2 at ±30 rad/s from Fig. 3
Fig. 17
Reference data, as well as the linear estimate of the x- and y-components of the angular velocity from the generated dataset with saturation limit 2 at ±30 rad/s from Fig. 3
Close modal
Fig. 18
Absolute error of the nonlinear estimate with respect to the reference angular velocities for the generated dataset with saturation limit 2 at ±30 rad/s from Fig. 3, thus two components at a time being saturated
Fig. 18
Absolute error of the nonlinear estimate with respect to the reference angular velocities for the generated dataset with saturation limit 2 at ±30 rad/s from Fig. 3, thus two components at a time being saturated
Close modal
Fig. 19
Comparison of different methods for orientation estimation for the generated dataset with two components at a time being saturated. The individual methods are defined in Table 2.
Fig. 19
Comparison of different methods for orientation estimation for the generated dataset with two components at a time being saturated. The individual methods are defined in Table 2.
Close modal
Fig. 20
Reference angular velocity as well as linear and nonlinear estimation of the z-component of the angular velocity from the field experiment with a rotating platform
Fig. 20
Reference angular velocity as well as linear and nonlinear estimation of the z-component of the angular velocity from the field experiment with a rotating platform
Close modal
Fig. 21
Error of the linear and nonlinear estimate with respect to the reference angular velocities for the field experiment with a rotating platform
Fig. 21
Error of the linear and nonlinear estimate with respect to the reference angular velocities for the field experiment with a rotating platform
Close modal
Fig. 22
Angular error of the linear and nonlinear estimate with respect to the reference angle
Fig. 22
Angular error of the linear and nonlinear estimate with respect to the reference angle
Close modal
Fig. 23
Comparison of different methods for orientation estimation for the field experiment with a rotating platform. The individual methods are defined in Table 2.
Fig. 23
Comparison of different methods for orientation estimation for the field experiment with a rotating platform. The individual methods are defined in Table 2.
Close modal
Fig. 24
Reference angular velocity as well as linear and nonlinear estimation of the z-component of the angular velocity from the acceleration phase of a snow avalanche
Fig. 24
Reference angular velocity as well as linear and nonlinear estimation of the z-component of the angular velocity from the acceleration phase of a snow avalanche
Close modal
Fig. 25
Error of the linear and nonlinear estimate with respect to the reference angular velocities for the acceleration phase of a snow avalanche
Fig. 25
Error of the linear and nonlinear estimate with respect to the reference angular velocities for the acceleration phase of a snow avalanche
Close modal
Fig. 26
Comparison of different methods for orientation estimation for the acceleration phase of a snow avalanche. The individual methods are defined in Table 2.
Fig. 26
Comparison of different methods for orientation estimation for the acceleration phase of a snow avalanche. The individual methods are defined in Table 2.
Close modal
Table 4

Mean and maximum errors of the orientation estimates for the simulated data experiment with two axes of the gyrometer in saturation shown in Fig. 19 

MethodMean error (rad)Max error (rad)
A4.205 03 × 10−22.037 95 × 10−1
B3.717 18 × 10−142.906 79 × 10−13
C2.906 04 × 10−19.694 08 × 10−1
D2.882 28 × 10−19.898 53 × 10−1
E8.456 75 × 10−13.019 77
MethodMean error (rad)Max error (rad)
A4.205 03 × 10−22.037 95 × 10−1
B3.717 18 × 10−142.906 79 × 10−13
C2.906 04 × 10−19.694 08 × 10−1
D2.882 28 × 10−19.898 53 × 10−1
E8.456 75 × 10−13.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 3.4×103 rad/s and 4.4×103 rad/s for the nonlinear and linear angular velocity estimate, respectively. Hence, the derived terminal angles are comparably accurate with an error of 1.6×102 rad and 2.1×102 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 103. 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.

Table 5

Mean and maximum errors of the orientation estimates of the rotating table experiment shown in Fig. 23 

MethodMean error (rad)Max error (rad)
A7.592 44 × 10−43.280 45 × 10−3
B2.239 75 × 10−41.855 92 × 10−3
C7.688 87 × 10−43.459 96 × 10−3
D2.408 54 × 10−41.437 39 × 10−3
E4.742 50 × 10−34.454 77 × 10−2
MethodMean error (rad)Max error (rad)
A7.592 44 × 10−43.280 45 × 10−3
B2.239 75 × 10−41.855 92 × 10−3
C7.688 87 × 10−43.459 96 × 10−3
D2.408 54 × 10−41.437 39 × 10−3
E4.742 50 × 10−34.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.

Table 6

Mean and maximum errors of the orientation estimates of the Nordkette snow avalanche shown in Fig. 26 

MethodMean error (rad)Max error (rad)
A0.124,760.420,10
B0.118,830.432,29
C0.112,270.376,26
D0.107,430.390,38
E1.170,362.910,06
MethodMean error (rad)Max error (rad)
A0.124,760.420,10
B0.118,830.432,29
C0.112,270.376,26
D0.107,430.390,38
E1.170,362.910,06

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.

Fig. 27
Saturated angular velocity as well as linear and nonlinear estimation of the z-component of the angular velocity from a snow avalanche with the gyrometer in saturation
Fig. 27
Saturated angular velocity as well as linear and nonlinear estimation of the z-component of the angular velocity from a snow avalanche with the gyrometer in saturation
Close modal

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×1012 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

2

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.

3

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

References

1.
Wall
,
J. H.
, and
Bevly
,
D. M.
et al
2005
, “
Characterization of Various IMU Error Sources and the Effect on Navigation Performance
,”
Proceedings of the 18th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS 2005)
, Long Beach, CA, Sept. 13–16, pp.
967
978
.
2.
Stančin
,
S.
, and
Tomažič
,
S.
,
2014
, “
Time-and Computation-Efficient Calibration of Mems 3d Accelerometers and Gyroscopes
,”
Sensors
,
14
(
8
), pp.
14885
14915
.10.3390/s140814885
3.
Renaudin
,
V.
,
Afzal
,
M. H.
, and
Lachapelle
,
G.
,
2010
, “
Complete Triaxis Magnetometer Calibration in the Magnetic Domain
,”
J. Sensors
,
2010
, pp.
1
10
.10.1155/2010/967245
4.
Kalman
,
R. E.
,
1960
, “
A New Approach to Linear Filtering and Prediction Problems
,”
ASME J. Basic Eng.
,
82
(
1
), pp.
35
45
.10.1115/1.3662552
5.
Madgwick
,
S.
,
2010
, “
An Efficient Orientation Filter for Inertial and Inertial/Magnetic Sensor Arrays
,”
Rep. x-io Univ. Bristol (UK)
,
25
, pp.
113
118
.
6.
Neurauter
,
R.
, and
Gerstmayr
,
J.
,
2023
, “
A Novel Motion-Reconstruction Method for Inertial Sensors With Constraints
,”
Multibody Syst. Dyn.
,
57
(
2
), pp.
181
209
.10.1007/s11044-022-09863-8
7.
Ojeda
,
L.
, and
Borenstein
,
J.
,
2007
, “
Personal Dead-Reckoning System for GPS-Denied Environments
,”
IEEE International Workshop on Safety, Security and Rescue Robotics
,
IEEE
, Rome, Italy, Sept. 27-29, pp.
1
6
.10.1109/SSRR.2007.4381271
8.
Zhang
,
R.
,
Yang
,
H.
,
Höflinger
,
F.
, and
Reindl
,
L. M.
,
2017
, “
Adaptive Zero Velocity Update Based on Velocity Classification for Pedestrian Tracking
,”
IEEE Sens. J.
,
17
(
7
), pp.
2137
2145
.10.1109/JSEN.2017.2665678
9.
Tan
,
C. H.
,
bin Shaiful
,
D. S.
,
Tang
,
E.
,
Khaw
,
J.-Y.
,
Soh
,
G. S.
, and
Foong
,
S.
,
2020
, “
Flydar: Magnetometer-Based High Angular Rate Estimation During Gyro Saturation for Slam
,”
IEEE International Conference on Robotics and Automation (ICRA)
,
IEEE
, Paris, France, May 31–Aug. 31, pp.
8532
8537
.10.1109/ICRA40945.2020.9197486
10.
Fischer
,
J.
, and
Rammer
,
L.
,
2010
, “
An Introduction to Inflow Avalanche Dynamics Measurements Using the Snowball Device
,”
16th International Snow Science Workshop (ISSW 2010)
,
Squaw Valley, CL
, Jan., pp.
738
741
.
11.
Vilajosana
,
I.
,
Llosa
,
J.
,
Schaefer
,
M.
,
Surin̂ach
,
E.
, and
Vilajosana
,
X.
,
2011
, “
Wireless Sensors as a Tool to Explore Avalanche Internal Dynamics: Experiments at the Weissflühjoch Snow Chute
,”
Cold Reg. Sci. Technol.
,
65
(
2
), pp.
242
250
.10.1016/j.coldregions.2010.09.011
12.
Dost
,
J. B.
,
Gronz
,
O.
,
Casper
,
M. C.
, and
Krein
,
A.
,
2020
, “
The Potential of Smartstone Probes in Landslide Experiments: How to Read Motion Data
,”
Nat. Hazards Earth Syst. Sci.
,
20
(
12
), pp.
3501
3519
.10.5194/nhess-20-3501-2020
13.
Kok
,
M.
, and
Schön
,
T. B.
,
2016
, “
Magnetometer Calibration Using Inertial Sensors
,”
IEEE Sens. J.
,
16
(
14
), pp.
5679
5689
.10.1109/JSEN.2016.2569160
14.
Winkler
,
R.
,
Fischer
,
J.-T.
,
Hergel
,
P.
,
Neuhauser
,
M.
,
Sovillia
,
B.
, and
Steinkogler
,
W.
,
2018
, “
Challenges and Limitations for in Situ Particle Tracking in Avalanches
,”
Proceedings of the International Snow Science Workshop (ISSW)
,
Innsbruck, Austria
, Squaw Valley, CF, Jan.
15.
Julier
,
S. J.
, and
Uhlmann
,
J. K.
,
1997
, “
New Extension of the Kalman Filter to Nonlinear Systems
,”
SPIE Proc.
, Vol.
3068
, pp.
182
193
.10.1117/12.280797
16.
Li
,
Q.
,
Li
,
R.
,
Ji
,
K.
, and
Dai
,
W.
,
2015
, “
Kalman Filter and Its Application
,”
8th International Conference on Intelligent Networks and Intelligent Systems (ICINIS)
, Tianjin, China, Nov., pp.
74
77
.
17.
Wahba
,
G.
,
1965
, “
A Least Squares Estimate of Satellite Attitude
,”
SIAM Rev.
,
7
(
3
), pp.
409
409
.10.1137/1007077
18.
Shuster
,
M. D.
, and
Oh
,
S. D.
,
1981
, “
Three-Axis Attitude Determination From Vector Observations
,”
J. Guid. Control
,
4
(
1
), pp.
70
77
.10.2514/3.19717
19.
Mortari
,
D.
,
1997
, “
ESOQ: A Closed-Form Solution to the Wahba Problem
,”
J. Astronaut. Sci.
,
45
(
2
), pp.
195
204
.10.1007/BF03546376
20.
Markley
,
F. L.
,
1993
, “
Attitude Determination Using Vector Observations: A Fast Optimal Matrix Algorithm
,”
J. Astronaut. Sci.
, 41(2), pp. 261–280.https://www.researchgate.net/publication/4718386_Attitude_Determination_from_Vector_Observations_A_Fast_Optimal_Matrix_Algorithm
21.
Mortari
,
D.
,
Markley
,
F. L.
, and
Singla
,
P.
,
2007
, “
Optimal Linear Attitude Estimator
,”
J. Guid., Control, Dyn.
,
30
(
6
), pp.
1619
1627
.10.2514/1.29568
22.
Madgwick
,
S. O.
,
Harrison
,
A. J.
, and
Vaidyanathan
,
R.
,
2011
, “
Estimation of IMU and MARG Orientation Using a Gradient Descent Algorithm
,”
IEEE International Conference on Rehabilitation Robotics
,
IEEE
, Zurich, Switzerland, June 29–July 1, pp.
1
7
.10.1109/ICORR.2011.5975346
23.
Potter
,
M. V.
,
Ojeda
,
L. V.
,
Perkins
,
N. C.
, and
Cain
,
S. M.
,
2019
, “
Effect of IMU Design on IMU-Derived Stride Metrics for Running
,”
Sensors
,
19
(
11
), p.
2601
.10.3390/s19112601
24.
Müller
,
A.
,
2010
, “
Approximation of Finite Rigid Body Motions From Velocity Fields
,”
ZAMM-J. Appl. Math. Mech./Z. Für Angew. Mathematik Und Mechanik: Appl. Math. Mech.
,
90
(
6
), pp.
514
521
.10.1002/zamm.200900383
25.
Müller
,
A.
,
2017
, “
Coordinate Mappings for Rigid Body Motions
,”
ASME J. Comput. Nonlinear Dyn.
,
12
(
2
), p.
021010
.10.1115/1.4034730
26.
Terze
,
Z.
,
Zlatar
,
D.
, and
Pandža
,
V.
,
2020
, “
Aircraft Attitude Reconstruction Via Novel Quaternion-Integration Procedure
,”
Aerosp. Sci. Technol.
,
97
, p.
105617
.10.1016/j.ast.2019.105617
27.
Murray
,
R. M.
,
Li
,
Z.
,
Sastry
,
S. S.
, and
Sastry
,
S. S.
,
1994
,
A Mathematical Introduction to Robotic Manipulation
,
CRC Press
, Boca Raton, FL.
28.
Holzinger
,
S.
, and
Gerstmayr
,
J.
,
2021
, “
Time Integration of Rigid Bodies Modelled With Three Rotation Parameters
,”
Multibody Syst. Dyn.
,
53
(
4
), pp.
345
378
.10.1007/s11044-021-09778-w
29.
Virtanen
,
P.
,
Gommers
,
R.
,
Oliphant
,
T. E.
,
Haberland
,
M.
,
Reddy
,
T.
,
Cournapeau
,
D.
,
Burovski
,
E.
, et al.,
2020
, “
SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python
,”
Nat. Methods
,
17
(
3
), pp.
261
272
.10.1038/s41592-019-0686-2
30.
Brüls
,
O.
, and
Cardona
,
A.
,
2010
, “
On the Use of Lie Group Time Integrators in Multibody Dynamics
,”
ASME J. Comput. Nonlinear Dyn.
,
5
(
3
), p.
031002
.10.1115/1.4001370
31.
Lourakis
,
M.
, and
Terzakis
,
G.
,
2018
, “
Efficient Absolute Orientation Revisited
,”
IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)
,
IEEE
, Madrid, Spain, Oct. 1–5, pp.
5813
5818
.10.1109/IROS.2018.8594296
32.
Neurauter
,
R.
,
Hergel
,
P.
, and
Gerstmayr
,
J.
,
2021
, “
Evaluation of Inertial Measurement Units for Short Time Motion Tracking
,”
ASME
Paper No. DETC2021-69604.10.1115/DETC2021-69604
33.
Terze
,
Z.
,
Müller
,
A.
, and
Zlatar
,
D.
,
2016
, “
Singularity-Free Time Integration of Rotational Quaternions Using Non-Redundant Ordinary Differential Equations
,”
Multibody Syst. Dyn.
,
38
(
3
), pp.
201
225
.10.1007/s11044-016-9518-7
34.
Holzinger
,
S.
,
Schöberl
,
J.
, and
Gerstmayr
,
J.
,
2020
, “
The Equations of Motion for a Rigid Body Using Non-Redundant Unified Local Velocity Coordinates
,”
Multibody Syst. Dyn.
,
48
(
3
), pp.
283
309
.10.1007/s11044-019-09700-5
35.
Arnold
,
M.
, and
Hante
,
S.
,
2017
, “
Implementation Details of a Generalized-α Differential-Algebraic Equation Lie Group Method
,”
ASME J. Comput. Nonlinear Dyn.
,
12
(
2
), p.
021002
.10.1115/1.4033441
36.
Gerstmayr
,
J.
,
2023
, “
Exudyn – A C++ based Python package for flexible multibody systems
,” Preprint, Research Square.
37.
Sonneville
,
V.
,
Cardona
,
A.
, and
Brüls
,
O.
,
2014
, “
Geometrically Exact Beam Finite Element Formulated on the Special Euclidean Group SE (3)
,”
Comput. Methods Appl. Mech. Eng.
,
268
, pp.
451
474
.10.1016/j.cma.2013.10.008