Numerical quaternion propagation and IMU attitude-state integration

Quaternion kinematics gives a continuous differential equation for attitude. A real navigation computer, flight controller, or simulation must convert that continuous equation into a discrete update executed at the IMU sample rate. This is the numerical quaternion propagation problem.

For an inertial measurement unit, the gyroscopes supply the angular-rate information used for the attitude prediction. The accelerometers do not normally determine the high-rate quaternion propagation directly. Instead, their measured specific force is rotated by the propagated attitude and used elsewhere in the navigation state equations, while accelerometer or gravity information may also be used in a later attitude-correction step.

This entry develops the exact finite-increment update, the small-angle approximation, forward Euler and Runge–Kutta propagation, gyro bias correction, delta-angle IMU data, coning effects, sign continuity, and the role of a quaternion inside a larger inertial-navigation or attitude-estimation state.

Convention declaration

This article continues the PhysicsLibrary convention used in Q11:

right-handed orthonormal frames;
Hamilton multiplication;
scalar-first quaternion display,

$\displaystyle q= \begin{bmatrix} q_w&q_x&q_y&q_z \end{bmatrix}^{T};$ (1)
the attitude quaternion is

$\displaystyle q\equiv{}^{I}q_B,$ (2)

mapping body-frame coordinates into inertial-frame coordinates;
the IMU gyroscope measurement is resolved in body axes unless stated otherwise.

With body-resolved angular velocity, Q11 derived

$\displaystyle \boxed{ \dot q=\frac12 q\,\omega^{B}, }$

(3)

where $\omega^{B}$ is embedded as a pure quaternion.

A sampled gyroscope model

A useful first-order gyroscope measurement model is

$\displaystyle \boxed{ \boldsymbol\omega_m^{B} = \boldsymbol\omega^{B} + \mathbf b_g + \mathbf n_g, }$

(4)

where $\mathbf b_g$ is gyro bias and $\mathbf n_g$ represents measurement noise. An estimator or navigation solution uses its current bias estimate $\widehat{\mathbf b}_g$ to form

$\displaystyle \boxed{ \boldsymbol\omega_c^{B} = \boldsymbol\omega_m^{B} - \widehat{\mathbf b}_g. }$

(5)

For a sample interval $\Delta t$ , the simplest delta-angle approximation is

$\displaystyle \boxed{ \Delta\boldsymbol\theta_k = \boldsymbol\omega_{c,k}^{B}\Delta t. }$

(6)

If the IMU itself reports integrated delta angles rather than angular rates, then the measured $\Delta\boldsymbol\theta_k$ may be used directly after applying the appropriate bias, scale-factor, and calibration corrections.

Exact quaternion increment for a constant sample rate

Assume that $\boldsymbol\omega_c^{B}$ is constant over one sample. Define

$\displaystyle \Delta\theta = \left\Vert \Delta\boldsymbol\theta \right\Vert,$

(7)

and, for nonzero $\Delta\theta$ ,

$\displaystyle \widehat{\mathbf u} = \frac{\Delta\boldsymbol\theta}{\Delta\theta}.$

(8)

The exact unit-quaternion increment is

$\displaystyle \boxed{ \delta q = \cos\frac{\Delta\theta}{2} + \widehat{\mathbf u} \sin\frac{\Delta\theta}{2}. }$

(9)

Equivalently, in scalar-first column form,

$\displaystyle \boxed{ \delta q = \begin{bmatrix} \cos(\Delta\theta/2)\\ [1mm] \dfrac{\sin(\Delta\theta/2)}{\Delta\theta} \Delta\boldsymbol\theta \end{bmatrix}. }$

(10)

Because the gyro rate is body-resolved, the increment multiplies the attitude state on the right:

$\displaystyle \boxed{ q_{k+1} = q_k\,\delta q_k. }$

(11)

This exponential update is exact when the angular velocity is constant over the sample interval.

Small-angle approximation

For IMU sampling at sufficiently high rate, the incremental rotation is usually small. Expanding the exact increment gives

$\displaystyle \cos\frac{\Delta\theta}{2} = 1-\frac{\Delta\theta^2}{8} +O(\Delta\theta^4),$

(12)

and

$\displaystyle \frac{\sin(\Delta\theta/2)}{\Delta\theta} = \frac12 - \frac{\Delta\theta^2}{48} + O(\Delta\theta^4).$

(13)

Thus

$\displaystyle \delta q = \begin{bmatrix} 1-\dfrac{\Delta\theta^2}{8}\\ [1mm] \l... ...\theta^2}{48} \right) \Delta\boldsymbol\theta \end{bmatrix}+ O(\Delta\theta^4),$

(14)

and the familiar first-order form is

$\displaystyle \boxed{ \delta q \approx \begin{bmatrix} 1\\ [1mm] \dfrac12\Delta\boldsymbol\theta \end{bmatrix}. }$

(15)

The first-order formula has two important applications, and they should not be confused.

Application 1: very small IMU propagation increments

For a sufficiently small gyro increment, one may use

$\displaystyle q_{k+1} \approx q_k \begin{bmatrix} 1\\ [1mm] \dfrac12\Delta\boldsymbol\theta_k \end{bmatrix},$

(16)

followed by normalization. This is computationally cheap, but the exact exponential increment is usually preferable when trigonometric evaluation is not a burden.

Application 2: local attitude errors

In an error-state estimator, the nominal attitude remains a full unit quaternion while a local attitude error is represented by the three-vector $\delta\boldsymbol\theta$ . From Q12,

$\displaystyle \boxed{ \delta q_e \approx \begin{bmatrix} 1\\ [1mm] \dfrac12\delta\boldsymbol\theta \end{bmatrix}. }$

(17)

This is a linearization of the error, not a replacement for the nominal quaternion state.

Figure summarizes the distinction between these two uses.

**Figure:** A typical IMU attitude-propagation flow. Body-frame gyro measurements are corrected for estimated bias, integrated to a delta angle, converted to an incremental quaternion, and multiplied on the right of the nominal attitude. The small-angle quaternion has two distinct uses: approximating a very small sample increment and representing a local three-component attitude error.

Forward Euler integration

One may discretize the continuous kinematic equation directly:

$\displaystyle q_{k+1}^{\mathrm{Euler}} = q_k + \frac12 q_k\omega_k^{B}\Delta t.$

(18)

This is precisely the first-order Taylor approximation to the exponential update. The difficulty is that the Euler step does not preserve unit norm exactly. Therefore implementations commonly apply

$\displaystyle \boxed{ q_{k+1} \leftarrow \frac{q_{k+1}^{\mathrm{Euler}}} {\left\Vert q_{k+1}^{\mathrm{Euler}}\right\Vert}. }$

(19)

Normalization corrects norm drift, but it does not recover all of the phase accuracy lost by a low-order integration method. At high dynamics or long sample intervals, exact or higher-order propagation is preferable.

Midpoint and exponential propagation

If the angular rate changes across the sample, a useful second-order strategy is to approximate the rate at the midpoint:

$\displaystyle \boldsymbol\omega_{k+1/2}^{B} \approx \frac12 \left( \boldsymbol\omega_k^{B} + \boldsymbol\omega_{k+1}^{B} \right).$

(20)

Then form

$\displaystyle \Delta\boldsymbol\theta_k = \boldsymbol\omega_{k+1/2}^{B}\Delta t$

(21)

and apply the exact quaternion increment of equation (10). This retains the unit-quaternion structure while improving the approximation of a varying rate.

Runge–Kutta propagation

The differential equation may also be integrated with a standard Runge–Kutta method. For example, fourth-order Runge–Kutta evaluates

$\displaystyle f(q,t) = \frac12q\,\omega^{B}(t)$

(22)

at four points within the interval. The resulting quaternion should normally be normalized after the update because generic numerical integration does not preserve the unit sphere exactly.

Runge–Kutta is valuable when angular rate is supplied by a continuous model. For sampled IMU data, however, an exponential or coning-corrected delta-angle update is often more natural because the sensor already provides discrete angular information.

A numerical IMU example

Suppose a body rotates at

$\displaystyle \boldsymbol\omega^{B} = \begin{bmatrix} 0&0&90 \end{bmatrix}^{T}$ deg/s

(23)

and the IMU is sampled at

$\displaystyle \Delta t=0.01$ s $\displaystyle .$

(24)

The delta angle per sample is

$\displaystyle \Delta\boldsymbol\theta = \begin{bmatrix} 0&0&0.9^\circ \end{bmatrix}^{T} = \begin{bmatrix} 0&0&0.015708 \end{bmatrix}^{T}$ rad

(25)

The exact increment is therefore approximately

$\displaystyle \delta q = \begin{bmatrix} 0.999969\ 0\ 0\ 0.007854 \end{bmatrix}.$

(26)

Applying this increment 100 times produces a $90^\circ$ rotation about the body $z$

axis, up to floating-point roundoff.

Bias error and attitude drift

Gyro bias integrates directly into attitude error. If an uncorrected constant bias $\delta\mathbf b_g$ is small, then over a short interval

$\displaystyle \delta\boldsymbol\theta \approx \delta\mathbf b_g\,t.$

(27)

Thus even a small gyro bias can create a steadily growing attitude error.

A typical estimator therefore includes gyro bias in the state:

$\displaystyle \boxed{ \mathbf x = \begin{bmatrix} q& \mathbf b_g^T& \cdots \end{bmatrix}^{T}. }$

(28)

A simple random-walk bias model is

$\displaystyle \dot{\mathbf b}_g = \mathbf n_{bg},$

(29)

or in discrete time,

$\displaystyle \mathbf b_{g,k+1} = \mathbf b_{g,k} + \mathbf w_{bg,k}.$

(30)

Small-angle error-state propagation

Consider a right, body-resolved attitude error as in Q12:

$\displaystyle q_{\mathrm{true}} = \widehat q\,\delta q_B.$

(31)

Let

$\displaystyle \delta q_B \approx \begin{bmatrix} 1\\ [1mm] \dfrac12\delta\boldsymbol\theta^{B} \end{bmatrix}.$

(32)

Define the gyro-bias estimation error as

$\displaystyle \delta\mathbf b_g = \mathbf b_g-\widehat{\mathbf b}_g.$

(33)

To first order, a common right-error model satisfies

$\displaystyle \boxed{ \delta\dot{\boldsymbol\theta}^{B} \approx - [\widehat{\bo... ...^{B}]_{\times} \delta\boldsymbol\theta^{B} - \delta\mathbf b_g - \mathbf n_g. }$

(34)

A first-order discrete approximation is

$\displaystyle \boxed{ \delta\boldsymbol\theta_{k+1}^{B} \approx \left( I_3 - [\... ...symbol\theta_k^{B} - \delta\mathbf b_{g,k}\Delta t - \mathbf n_{g,k}\Delta t. }$

(35)

The signs and multiplication side are convention dependent. The state, injection rule, and covariance model must all use the same error convention.

Connection to the rest of an IMU navigation state

In strapdown navigation, the attitude update is only one part of the state propagation. The accelerometer measures specific force in body axes, $\mathbf f^B$ . The propagated quaternion supplies the rotation

$\displaystyle \mathbf f^{I} = R(q)\mathbf f^{B},$

(36)

which then enters the velocity and position equations. In a simple inertial frame model,

$\displaystyle \dot{\mathbf v}^{I} = \mathbf f^{I} + \mathbf g^{I}.$

(37)

More complete navigation equations add Earth rotation, transport rate, gravity models, and sensor calibration terms. The important structural point is that gyro integration propagates the quaternion used to orient the accelerometer measurement.

Multiple increments and coning

Rotations do not commute. Therefore several IMU delta-angle samples cannot in general be replaced by a simple vector sum.

For two small body-frame increments applied in chronological order,

$\displaystyle q_{k+2} = q_k\, \delta q_1 \delta q_2.$

(38)

Using the first-order quaternion form for each increment and retaining the leading cross term gives the equivalent small rotation vector

$\displaystyle \boxed{ \Delta\boldsymbol\theta_{\mathrm{eq}} \approx \Delta\bold... ...theta_2 + \frac12 \Delta\boldsymbol\theta_1 \times \Delta\boldsymbol\theta_2. }$

(39)

The cross-product term is the first indication of the coning effect. Practical strapdown algorithms use more accurate coning corrections that depend on how the gyro samples or delta angles are accumulated within the update interval.

Sign continuity

The quaternions $q$ and $-q$ represent the same physical orientation. A numerical time history may therefore contain an unnecessary sign jump if a conversion algorithm chooses the opposite representative at one sample.

A common continuity rule is

$\displaystyle \boxed{ \text{if } q_k^T q_{k-1}<0, \qquad q_k\leftarrow -q_k. }$

(40)

This does not change the physical attitude. It simply chooses the quaternion representative closest to the previous sample, which is useful for plotting, interpolation, and finite differencing.

Choosing a propagation method

For many IMU applications the following hierarchy is useful:

use the exact exponential increment when one rate or delta-angle vector is representative of the sample interval;
use a midpoint or coning-corrected increment when angular rate changes significantly within the interval;
use Runge–Kutta when integrating a continuous angular-rate model or when the system dynamics already use a compatible higher-order integrator;
reserve the first-order small-angle update mainly for very small increments, linearized error-state equations, and analytical derivations.

Common pitfalls

The most common implementation errors are:

multiplying a body-resolved increment on the left instead of the right;
applying the small-angle approximation to a large sample rotation;
subtracting no gyro bias, or using the wrong sign for the bias estimate;
assuming normalization makes a low-order integrator exact;
summing finite delta-angle vectors while ignoring noncommutativity and coning;
confusing the small three-component attitude error with the full nominal quaternion state;
rotating accelerometer specific force with a stale or oppositely defined attitude quaternion;
treating the sign change $q\leftrightarrow -q$ as a physical attitude jump.

What comes next

Q14 consolidates the quaternion sequence into a larger graded exercise bank. The numerical ideas developed here also provide the direct bridge into strapdown inertial navigation, multiplicative attitude filters, and quaternion-based simulation states.

Sources and engineering context

The continuous kinematic equation follows from Q11 and the small-error representation follows from Q12. The discrete exponential update is the finite-rotation form of that same kinematics. Gyro bias, delta-angle integration, coning, and attitude-state propagation are standard concerns in strapdown inertial navigation and spacecraft attitude computation.

License

This article is an original synthesis prepared for PhysicsLibrary and intended for release under CC BY-SA 4.0.

"numerical quaternion propagation and IMU attitude state integration" is owned by bloftin.
(view preamble)