M6-T4 — INS Error Dynamics and the Schuler Oscillation

Learning objectives

After this lesson, you should be able to:

• Write the error equations of one horizontal channel and the vertical channel.
• Derive the Schuler frequency and the position response to an accelerometer and a gyro bias.
• Derive the exponential growth of the vertical channel and its time constant.
• State the time scales on which the Schuler effect does and does not matter.

The intuition

UAV M2-T8: Dead reckoning and inertial drift states that inertial position error grows with time. This lesson shows the shape of the growth and its causes, for one channel, with equations small enough to solve in closed form and check against a simulation.

Why the horizontal error does not grow without bound

Suppose the computer thinks it has moved a little farther east than it has. It then also thinks it has travelled over a curved Earth, so the local vertical it computes has turned by the position error divided by the radius of the Earth. The platform is tilted by that angle relative to the true local level. A tilt makes the gravity vector leak into the horizontal accelerometer, and the leakage is a restoring force: it pulls the computed position back. A mass on a spring is the picture, and the spring constant is $g/R$. The frequency is $\omega_s = \sqrt{g/R} = 1.24\times10^{-3}$ rad/s, a period of 84.4 minutes, the period of a pendulum whose length equals the radius of the Earth. This is the Schuler oscillation.

What the sensor errors do to it

An accelerometer bias is a constant force on that spring, so the error oscillates about a displaced centre with amplitude $b/\omega_s^2$: 0.1 mg gives a peak error of 1.27 km (at 42 min) and a return to zero at 84 min. A gyro bias is a constant tilt rate: it makes the error grow steadily at $\varepsilon R$ on average (1.1 km per hour for 0.01°/h) with the oscillation riding on it. In both cases, over times much shorter than the period, the coupling has not yet acted, and the errors follow the polynomial growth of the next lesson: $\tfrac{1}{2}bt^2$ for the accelerometer, $\tfrac{1}{6}g\varepsilon t^3$ for the gyro.

The vertical channel has no restoring force

Gravity falls with height, so a computed height that is too high gives a computed gravity that is too small, a smaller restoring pull than needed and a net upward acceleration: positive feedback. The vertical channel is unstable, with a growth rate $\sqrt{2g/R} = 1.76\times10^{-3}$ per second: an e-folding every 9.5 minutes, a doubling every 6.6. A 0.1 mg bias gives 1.8 m after 60 s, 194 m after 10 min and 88 km after an hour. Every inertial system needs a barometer or a radar altimeter to hold the vertical channel.

The Schuler coupling makes horizontal inertial errors oscillate with an 84-minute period set by the Earth's radius alone. The vertical channel has positive feedback and diverges with a time constant of 9.5 minutes. On a 60 s hop neither matters, and the errors are the simple polynomials.

Definitions

Schuler frequency $\omega_s$

$\omega_s = \sqrt{g/R}$ = 1.24×10⁻³ rad/s, period $T_s = 2\pi\sqrt{R/g}$ = 84.4 min.

Position error $\delta x$, velocity error $\delta v$, tilt $\psi$

The errors of one horizontal channel. The tilt is the angle between the computed and the true local vertical.

Accelerometer bias $b_a$, gyro bias $\varepsilon$

The two sensor errors that drive the channel.

Vertical channel

The height and vertical velocity errors; coupled to gravity variation with height, not to a tilt.

Derivation: the Schuler frequency, and the responses

Eliminate the tilt from the horizontal channel.

1. $$\ddot{\delta x} = -g\,\psi + b_a$$

   Differentiate the position error and use the velocity equation.

2. $$\dddot{\delta x} = -g\,\dot{\psi} = -g\left(\dfrac{\delta v}{R} + \varepsilon\right)$$

   Differentiate again and substitute the tilt rate. The accelerometer bias is constant and drops out.

3. $$\ddot{\delta v} + \dfrac{g}{R}\,\delta v = -g\,\varepsilon$$

   Write the equation for $\delta v = \dot{\delta x}$. It is a harmonic oscillator with frequency $\omega_s^2 = g/R$, driven by the gyro bias.

4. $$\omega_s = \sqrt{g/R} = \sqrt{9.81/6.371\times10^6} = 1.241\times10^{-3}\ \mathrm{rad/s}$$

   The frequency of the free oscillation. Period $2\pi/\omega_s$ = 5063 s = 84.4 min, whatever the vehicle.

The responses to the two biases

1. $$\delta x(t) = \dfrac{b_a}{\omega_s^2}\bigl(1 - \cos\omega_s t\bigr)$$

   The solution for an accelerometer bias from zero initial errors. For $\omega_s t \ll 1$ it becomes $\tfrac12 b_a t^2$: the coupling has not acted.

2. $$\delta x(t) = \varepsilon R\left(t - \dfrac{\sin\omega_s t}{\omega_s}\right)$$

   The solution for a gyro bias: a steady drift at $\varepsilon R$ with an oscillation of amplitude $\varepsilon R/\omega_s$. For short times it becomes $\tfrac16 g\varepsilon t^3$.

3. $$\delta h(t) \propto e^{\lambda t}, \quad \lambda = \sqrt{2g/R} = 1.755\times10^{-3}\ \mathrm{s^{-1}}$$

   The vertical channel has the characteristic equation $s^2 = 2g/R$, with a positive real root: e-folding time 9.5 min, doubling time 6.6 min.

The coupling between position error, tilt and gravity turns the horizontal channel into an undamped oscillator with the Schuler frequency, and the vertical channel into an unstable first-order growth. Both come from the same g/R, with different signs.

A worked example

One channel simulated with the three-equation error model and compared with the closed forms.

1. Frequencies. The Schuler frequency is 1.241×10⁻³ rad/s, a period of 84.4 min. The vertical growth rate is 1.755×10⁻³ s⁻¹, an e-folding of 9.5 min and a doubling of 6.6 min.
2. Accelerometer bias of 0.1 mg. Position error 1.76 m at 60 s (simulated and closed form agree, the parabola gives 1.77), 168.6 m at 600 s (the parabola gives 176.6), 1274 m at half a period (the parabola 3140), and back to 0.01 m at one period.
3. Gyro bias of 0.01°/h. Position error 0.02 m at 60 s, 16.7 m at 600 s, 1353 m at one hour and 1564 m at 5060 s. The long-time drift is $\varepsilon R$ = 1112 m per hour: a nautical mile per hour needs 0.0167°/h.
4. Vertical, 0.1 mg bias. Height error 1.8 m at 60 s (the parabola 1.8), 193 m at 600 s (177), 3438 m at 1800 s (1589) and 88 km at 3600 s (6.4 km): the exponent takes over after about ten minutes.
5. Reading for Hopper. A hop lasts about a minute, 1/84 of a Schuler period. Over that time the errors are the parabola and the cubic to 0.1%, and the Schuler and vertical-channel instabilities have not started.

import math
import gnckit as g

Re, G0 = 6371000.0, 9.81
ws = math.sqrt(G0 / Re)
print(f'Schuler frequency sqrt(g/R) = {ws:.3e} rad/s, period {2 * math.pi / ws / 60:.1f} min')
wv = math.sqrt(2 * G0 / Re)
print(f'vertical channel: growth rate sqrt(2g/R) = {wv:.3e} per s, e-folding {1 / wv / 60:.1f} min, doubling {math.log(2) / wv / 60:.1f} min')

# one horizontal channel error equations: d(dx)' = dv, dv' = -g psi + ba, psi' = dv / R + eps
def run(ba, eps, T, dt=0.5):
    def f(t, s):
        dx, dv, psi = s
        return [dv, -G0 * psi + ba, dv / Re + eps]
    return g.integrate(f, [0.0, 0.0, 0.0], T, dt)

mg = 1e-3 * G0
ba = 0.1 * mg
print('accelerometer bias 0.1 mg: position error over time')
for T in (60.0, 600.0, 2530.0, 5060.0):
    ts, xs = run(ba, 0.0, T)
    an = ba / ws ** 2 * (1 - math.cos(ws * T))
    print(f'  t = {T:6.0f} s: simulated {xs[-1][0]:10.2f} m, formula (b/w^2)(1 - cos wt) = {an:10.2f} m, early-time 0.5 b t^2 = {0.5 * ba * T * T:10.2f} m')
eps = 0.01 * math.pi / 180 / 3600
print('gyro bias 0.01 deg/h:')
for T in (60.0, 600.0, 3600.0, 5060.0):
    ts, xs = run(0.0, eps, T)
    an = eps * Re * (T - math.sin(ws * T) / ws)
    print(f'  t = {T:6.0f} s: simulated {xs[-1][0]:10.2f} m, formula eps R (t - sin(wt)/w) = {an:10.2f} m, early-time g eps t^3 / 6 = {G0 * eps * T ** 3 / 6:10.2f} m')
print(f'  the long-time drift is eps R = {eps * Re * 3600:.1f} m per hour: a nautical mile per hour corresponds to {1852 / 3600 / Re / (math.pi / 180) * 3600:.4f} deg/h')

# vertical channel with a 0.1 mg accelerometer bias
def runv(ba, T, dt=0.5):
    def f(t, s):
        return [s[1], 2 * G0 / Re * s[0] + ba]
    return g.integrate(f, [0.0, 0.0], T, dt)
for T in (60.0, 600.0, 1800.0, 3600.0):
    ts, xs = runv(ba, T)
    print(f'vertical channel, accelerometer bias 0.1 mg, t = {T:5.0f} s: height error {xs[-1][0]:10.1f} m (0.5 b t^2 = {0.5 * ba * T * T:9.1f} m)')

Schuler frequency sqrt(g/R) = 1.241e-03 rad/s, period 84.4 min
vertical channel: growth rate sqrt(2g/R) = 1.755e-03 per s, e-folding 9.5 min, doubling 6.6 min
accelerometer bias 0.1 mg: position error over time
  t =     60 s: simulated       1.76 m, formula (b/w^2)(1 - cos wt) =       1.76 m, early-time 0.5 b t^2 =       1.77 m
  t =    600 s: simulated     168.57 m, formula (b/w^2)(1 - cos wt) =     168.57 m, early-time 0.5 b t^2 =     176.58 m
  t =   2530 s: simulated    1274.20 m, formula (b/w^2)(1 - cos wt) =    1274.20 m, early-time 0.5 b t^2 =    3139.64 m
  t =   5060 s: simulated       0.01 m, formula (b/w^2)(1 - cos wt) =       0.01 m, early-time 0.5 b t^2 =   12558.57 m
gyro bias 0.01 deg/h:
  t =     60 s: simulated      -0.02 m, formula eps R (t - sin(wt)/w) =       0.02 m, early-time g eps t^3 / 6 =       0.02 m
  t =    600 s: simulated     -16.65 m, formula eps R (t - sin(wt)/w) =      16.65 m, early-time g eps t^3 / 6 =      17.12 m
  t =   3600 s: simulated   -1353.42 m, formula eps R (t - sin(wt)/w) =    1353.42 m, early-time g eps t^3 / 6 =    3698.28 m
  t =   5060 s: simulated   -1563.98 m, formula eps R (t - sin(wt)/w) =    1563.98 m, early-time g eps t^3 / 6 =   10269.38 m
  the long-time drift is eps R = 1111.9 m per hour: a nautical mile per hour corresponds to 0.0167 deg/h
vertical channel, accelerometer bias 0.1 mg, t =    60 s: height error        1.8 m (0.5 b t^2 =       1.8 m)
vertical channel, accelerometer bias 0.1 mg, t =   600 s: height error      193.5 m (0.5 b t^2 =     176.6 m)
vertical channel, accelerometer bias 0.1 mg, t =  1800 s: height error     3437.8 m (0.5 b t^2 =    1589.2 m)
vertical channel, accelerometer bias 0.1 mg, t =  3600 s: height error    87953.5 m (0.5 b t^2 =    6356.9 m)

In robotics

The Schuler tuning and the vertical instability decide how an INS is built and aided.

• Long flights. For aircraft and ships the 84-minute oscillation is the dominant feature of the free-inertial error: the position error of a good gyro is a drift of $\varepsilon R$ with a Schuler ripple, and the typical specification is "1 nautical mile per hour" (0.0167°/h).
• Short flights. For a hop, a launch or a drone mission of minutes the Schuler effect does not act, and the polynomial laws (next lesson) are what size the sensors.
• Vertical aiding. A barometric or radar altimeter is always blended with the vertical channel. The simplest form is a third-order feedback loop whose gains are chosen to damp the exponential, as used in classical INS.
• Alignment error. An initial tilt of $\psi_0$ starts the oscillation with an amplitude $g\psi_0/\omega_s^2$ of the same form, which is why alignment accuracy matters as much as sensor bias.
• Aided systems. When a position fix is available the filter of module 7 estimates the tilt and bias states and cancels the oscillation. The Schuler dynamics are part of the model the filter uses.

The lab M6-L2: Mechanize and Watch It Drift reproduces these curves from a planar mechanization. The next two lessons align the system and budget the errors.

Common misconceptions

Key takeaways

1. One horizontal channel: $\dot{\delta x} = \delta v$, $\dot{\delta v} = -g\psi + b_a$, $\dot\psi = \delta v/R + \varepsilon$.
2. It oscillates at $\omega_s = \sqrt{g/R}$: period 84.4 min for any vehicle.
3. An accelerometer bias gives $\delta x = (b_a/\omega_s^2)(1 - \cos\omega_s t)$ (peak 1.27 km for 0.1 mg); a gyro bias gives a drift $\varepsilon R$ (1.1 km/h for 0.01°/h).
4. For times much shorter than the period the errors reduce to $\tfrac12 b_a t^2$ and $\tfrac16 g\varepsilon t^3$.
5. The vertical channel is unstable: $\ddot{\delta h} = (2g/R)\delta h + b_a$, e-folding 9.5 min; it needs an altimeter.
6. On a 60 s hop none of the Schuler and vertical effects have started.

References and further reading

Each source is marked free or paid. Free to read does not mean free to reuse: the course resources page lists every source with its licence.

• P. D. Groves, Principles of GNSS, Inertial, and Multisensor Integrated Navigation Systems, 2nd edition, Artech House, 2013, chapter 5 on inertial error propagation. Paid.
• D. H. Titterton and J. L. Weston, Strapdown Inertial Navigation Technology, 2nd edition, IET, 2004, the error-model chapters. Paid.
• M. Schuler, "Die Störung von Pendel- und Kreiselapparaten durch die Beschleunigung des Fahrzeuges", Physikalische Zeitschrift 24, 1923, the original paper (in German). Free to read in scanned journal archives.
• MIT 16.30 Feedback Control Systems (OpenCourseWare) , for the oscillator and feedback background used here. Free; licence CC BY-NC-SA 4.0.

Knowledge check

Three questions. Write an answer, then reveal the explanation. Your answers stay in your browser.

What’s next

Related content

Nothing below is required to finish this lesson.

• Before
  M6-T3: Strapdown Mechanization

  Write the attitude, velocity and position updates in the local-level frame, including gravity, Earth rate and transport rate, and say when coning and sculling corrections matter.

• Lab
  M6-L2: Mechanize and Watch It Drift

  Write a planar strapdown mechanization, inject a known accelerometer bias and gyro bias, and match the growth of the position error to the error equations.

Other courses on this site treat neighbouring topics from their own angle:

• UAV M2-T7: Inertial Navigation Systems
• UAV M2-T8: Dead reckoning and inertial drift
• UAV M3-T7: IMU noise and bias modelling
• UAV M3-T8: IMU preintegration
• UAV M1-T8: IMU, barometer, magnetometer and GNSS
• Probability M2-T4: Distributions of Real Sensor Noise
• Probability M5-T4: Identifying Noise Parameters from Data
• Linear Algebra M10-T6: Small-Angle Approximations

This course leaves the following to the pages that already teach them:

• What an INS is and why it drifts, without the equations
• IMU preintegration for factor-graph estimators
• Wheel odometry and its motion models