INS error propagation — error state equations: Exercises

Prerequisites you may want open: Strapdown INS Mechanization Equations, Gyro and Accelerometer Error Models, Direction Cosine Matrix and Small-Angle Rotations, Kalman Filter for INS-GPS Integration, Schuler Tuning and Oscillation.

Definition Notation you need before starting

Five pieces of shorthand recur below. Learn them once here.

1. The NED frame and its "Down" sign. We work in the local North-East-Down frame: the first axis points North, the second East, and the third points straight down toward Earth's centre. So a vector that points upward (away from the ground) has a negative third component. Gravity, which pulls down, has a positive Down component ( $+ g$ ). The accelerometers, however, feel the support force holding you up, so the specific force $f^{n}$ points up — its Down component is $- g$ . Whenever you see " $f_{D}^{n} = - g$ ", read it as "the specific force points up with size $g$ ."

2. The skew-symmetric bracket $[a \times]$ . For any vector $a = (a_{x}, a_{y}, a_{z})$ , the symbol $[a \times]$ is the $3 \times 3$ matrix that does the cross product by matrix multiplication: $[a \times] b = a \times b$ for every $b$ . Explicitly $[a \times] = 0 a_{z} - a_{y} - a_{z} 0 a_{x} a_{y} - a_{x} 0 .$ Why we use it: the error equations are linear in $δ x$ , and "cross with $a$ " is a linear operation, so packaging it as a matrix lets us drop every coupling neatly into the matrix $F$ . Reading $- [f^{n} \times] ψ$ just means " $- f^{n} \times ψ$ ."

3. The two "how the world turns" vectors $ω_{i e}$ and $ω_{e n}$ . These are small angular-rate vectors (units rad/s) that keep track of the fact that the navigation frame is not an inertial frame — it turns as the Earth spins and as the vehicle travels over the curved Earth. Read the subscripts as "from-frame to frame, resolved in-frame": $i$ = the fixed inertial frame (non-rotating, star-fixed), $e$ = the Earth frame (rotates with the planet), $n$ = the local nav (NED) frame.

$ω_{i e}$ = rate of the Earth frame relative to inertial — i.e. Earth's spin, magnitude $ω_{i e} = 7.292 \times 1 0^{- 5}$ rad/s (one turn per day).
$ω_{e n}$ = rate of the nav frame relative to the Earth — the "transport rate", how fast your local North-East-Down triad swings around as you fly across the globe (bigger when you move fast). It is $0$ if you sit still.
Their sum $ω_{in} = ω_{i e} + ω_{e n}$ is the total turn rate of the nav frame relative to inertial. This is the vector that carries the tilt vector $ψ$ along in the attitude equation.

4. The rotation matrix $C_{b}^{n}$ . The INS senses in the body frame (the axes bolted to the vehicle, superscript $b$ ); we navigate in the nav frame (superscript $n$ ). $C_{b}^{n}$ is the $3 \times 3$ direction-cosine matrix (see Direction Cosine Matrix and Small-Angle Rotations) that rotates a vector from body to nav coordinates: $v^{n} = C_{b}^{n} v^{b}$ . It is a pure rotation (orthogonal, determinant $+ 1$ ), so it only re-orients a vector — never stretches it. That is why sensor errors measured in body axes ( $δ f^{b}, δ ω^{b}$ ) get multiplied by $C_{b}^{n}$ before they can enter nav-frame velocity or attitude equations.

5. The input pair $w$ and $G$ , with block sizes. The sensor errors are collected into the 6×1 input vector $w = [δ f^{b}, δ ω^{b}]^{T}$ (a 3-vector accelerometer error stacked on a 3-vector gyro error, both in body frame). The state $δ x = [δ r, δ v, ψ]^{T}$ is 9×1 (three 3-vectors). So in $\dot{δ x} = F δ x + G w$ , $F$ is $9 \times 9$ and $G$ is $9 \times 6$ . Writing $G$ in $3 \times 3$ blocks, $G = 0_{3 \times 3} C_{b}^{n} 0_{3 \times 3} 0_{3 \times 3} 0_{3 \times 3} - C_{b}^{n},$ which is the compact " $diag (0, C_{b}^{n}, - C_{b}^{n})$ " — accelerometer error (left column) feeds the velocity rows, gyro error (right column) feeds the attitude rows, nothing feeds position. Sensor errors live in $w$ — never inside the state $δ x$ .

Level 1 — Recognition

Here we only name things. No algebra — just read the model and point to the right piece.

Recall Solution L1.1

What: $ψ$ is the attitude/tilt error, measured in radians. Picture: the INS carries a computed North-East-Down frame. If its bookkeeping is slightly wrong, that computed frame is rotated away from the true frame by a tiny turn. $ψ$ is the little arrow (rotation vector) describing that turn — its length is the tilt angle, its direction is the axis you'd twist about. $δ r$ is in metres, $δ v$ is in metres/second, $ψ$ is in radians.

Recall Solution L1.2

Position ← velocity: the identity block $I$ in the top row, second column. It encodes $\dot{δ r} = δ v$ exactly. Attitude → velocity: the $- [f^{n} \times]$ block (middle row, third column). A tilt $ψ$ mis-resolves the big specific-force vector $f^{n}$ , and that mis-projection appears as a velocity error. Why read $F$ this way: row $=$ "who am I", column $=$ "who pushes me". Entry $F_{ij}$ tells you how strongly state $j$ drives the rate of state $i$ . (The $- [(2 ω_{i e} + ω_{e n}) \times]$ block on the velocity diagonal is the Coriolis/transport reshuffling of velocity error as the frame turns.)

Recall Solution L1.3

Meaning of the symbols: $ω_{in} = ω_{i e} + ω_{e n}$ is the total turn rate of the nav frame relative to inertial (Earth spin plus transport rate); $C_{b}^{n}$ rotates the body-frame gyro error into nav coordinates.

Transport-rate coupling $- ω_{in} \times ψ$ : as the vehicle moves over the rotating Earth, the nav frame itself turns at rate $ω_{in}$ , and the tilt vector is carried along by that turn.
Gyro error $- C_{b}^{n} δ ω_{ib}^{b}$ : the gyros' drift and noise, rotated from body frame into nav frame. This is the fresh injection of error every instant — the root cause of long-term drift.

Level 2 — Application

Now we integrate. Every answer here is "count the integrations, get the power of $t$ ."

Recall Solution L2.1

Step 1 — convert the bias. $b_{a} = 300 \times 9.8 \times 1 0^{- 6} = 2.94 \times 1 0^{- 3} m/s^{2}$ . Step 2 — velocity (first integral). $\dot{δ v} = b_{a}$ constant $\Rightarrow δ v (t) = b_{a} t$ . $δ v (120) = 2.94 \times 1 0^{- 3} \times 120 = 0.3528 m/s$ . Why linear? One integration of a constant gives a straight ramp in $t$ . Step 3 — position (second integral). $δ r (t) = \frac{1}{2} b_{a} t^{2}$ . $δ r (120) = \frac{1}{2} \times 2.94 \times 1 0^{- 3} \times 12 0^{2} = 21.168 m$ . Why quadratic? Position is the second integral of the source, so $t^{2}$ .

Recall Solution L2.2

Step 1 — convert drift to rad/s. $ϵ = 0.01 \times \frac{π}{180} = 1.7453 \times 1 0^{- 4} rad/s$ . Step 2 — tilt (first integral of drift). $ψ (t) = ϵ t$ , so $ψ (60) = 1.7453 \times 1 0^{- 4} \times 60 = 1.0472 \times 1 0^{- 2}$ rad ( $\approx 0. 6^{\circ}$ ). Step 3 — velocity (second integral overall). $\dot{δ v} = - g ϵ t \Rightarrow δ v (t) = - \frac{1}{2} g ϵ t^{2}$ . $δ v (60) = - \frac{1}{2} \times 9.8 \times 1.7453 \times 1 0^{- 4} \times 6 0^{2} = - 3.0788 m/s$ . Step 4 — position (third integral overall). $δ r (t) = - \frac{1}{6} g ϵ t^{3}$ . $δ r (60) = - \frac{1}{6} \times 9.8 \times 1.7453 \times 1 0^{- 4} \times 6 0^{3} = - 61.576 m$ . Why cubic? drift → tilt (1) → velocity (2) → position (3): three integrations, hence $t^{3}$ .

Recall Solution L2.3

Accel channel: $b_{a} = 100 \times 9.8 \times 1 0^{- 6} = 9.8 \times 1 0^{- 4} m/s^{2}$ . $δ r_{a} = \frac{1}{2} b_{a} t^{2} = \frac{1}{2} \times 9.8 \times 1 0^{- 4} \times 10 0^{2} = 4.9 m$ . Gyro channel: $ϵ = 1.7453 \times 1 0^{- 4} rad/s$ . $δ r_{g} = \frac{1}{6} g ϵ t^{3} = \frac{1}{6} \times 9.8 \times 1.7453 \times 1 0^{- 4} \times 10 0^{3} = 285.06 m$ . Ratio: $285.06/4.9 \approx 58.2$ . The gyro channel dominates by roughly a factor of $58$ . Why: the extra integration ( $t^{3}$ vs $t^{2}$ ) plus the gravity gain $g$ make gyro drift the long-run killer.

Level 3 — Analysis

Now we reason about the couplings and signs, not just plug numbers.

Recall Solution L3.1

The cross product $f^{n} \times ψ$ with $f^{n} = (0, 0, - g)$ and $ψ = (0, ψ_{E}, 0)$ : $f^{n} \times ψ = i 00 j 0 ψ_{E} k - g 0 = i (0 \cdot 0 - (- g) ψ_{E}) - j (0) + k (0) = (g ψ_{E}, 0, 0) .$ So $- [f^{n} \times] ψ = - (g ψ_{E}, 0, 0)$ : it drives the North velocity error, $\dot{δ v}_{N} = - g ψ_{E}$ . Picture / why: an East-axis tilt "leans" the wrongly-oriented North accelerometer into the gravity vector, so a slice of $g$ leaks into North acceleration. See the tilt figure below.

Recall Solution L3.2

Compute the coupling angle over the flight: $ω_{i e} t = 7.292 \times 1 0^{- 5} \times 120 = 8.75 \times 1 0^{- 3}$ rad $\approx 0. 5^{\circ}$ . Because this is $≪ 1$ , the Earth-rate rotation reshuffles the error vector by under half a degree over the whole flight — a $< 1%$ effect on the dominant channels. Conclusion: yes, drop $ω_{i e}$ (and the comparably small $ω_{e n}$ ) for a short flight; keep only the $- [f^{n} \times]$ coupling and the sensor-error inputs. This is the 80/20 of short-flight INS error.

Recall Solution L3.3

Close the loop. Differentiate $\dot{δ v} = - g ψ$ once more and substitute $\dot{ψ} = \frac{1}{R} δ v$ : $\ddot{δ v} = - g \dot{ψ} = - g \cdot \frac{1}{R} δ v = - \frac{g}{R} δ v .$ This is the equation of simple harmonic motion, $\ddot{δ v} + ω_{s}^{2} δ v = 0$ with $ω_{s} = g / R$ . Since $δ r$ is the integral of $δ v$ , position error oscillates at the same frequency. Why oscillation, not a ramp: gravity closes the loop with a negative sign (the $- g$ in $\dot{δ v}$ combined with the $+ \frac{1}{R}$ feedback), exactly like a spring restoring force $- k x$ . Without this feedback ( $\frac{1}{R} \to 0$ ) the loop is open and the error ramps as $t^{3}$ ; with it, the error swings back — the Schuler oscillation at $\approx 84$ min. A constant accel bias still adds a steady offset the loop cannot null, so the real error is a bounded swing plus a slow ramp.

Level 4 — Synthesis

Assemble whole pieces of the model and run them.

Recall Solution L4.1

Reduced system: $\dot{δ r} = δ v, \dot{δ v} = - g ψ, \dot{ψ} = ϵ .$ Integrate top-down (each from rest):

$ψ (t) = ϵ t$ — first integral of the constant drift.
$\dot{δ v} = - g ϵ t \Rightarrow δ v (t) = - \frac{1}{2} g ϵ t^{2}$ — second integral overall.
$\dot{δ r} = - \frac{1}{2} g ϵ t^{2} \Rightarrow δ r (t) = - \frac{1}{6} g ϵ t^{3}$ — third integral overall. Result: $δ r (t) = - \frac{1}{6} g ϵ t^{3}$ — the cubic law from first principles, matching Example 2 of the parent note. Three integrations (drift→tilt→velocity→position) give the $t^{3}$ .

Recall Solution L4.2

Symbolic: $ψ = ϵ t$ , so $\dot{δ v} = b_{a} - g ϵ t \Rightarrow δ v = b_{a} t - \frac{1}{2} g ϵ t^{2}$ , and $δ r (t) = \frac{1}{2} b_{a} t^{2} - \frac{1}{6} g ϵ t^{3} .$ Numbers: $b_{a} = 200 \times 9.8 \times 1 0^{- 6} = 1.96 \times 1 0^{- 3} m/s^{2}$ ; $ϵ = 0.005 \times π /180 = 8.7266 \times 1 0^{- 5} rad/s$ .

Accel part: $\frac{1}{2} \times 1.96 \times 1 0^{- 3} \times 9 0^{2} = 7.938 m$ .
Gyro part: $\frac{1}{6} \times 9.8 \times 8.7266 \times 1 0^{- 5} \times 9 0^{3} = 103.90 m$ .
$δ r (90) = 7.938 - 103.90 = - 95.96 m$ . Read it: even a smaller drift dominates by $t = 90$ s because of the cubic power — synthesis confirms the gyro is the binding spec.

Level 5 — Mastery

Limits, degenerate cases, and design reasoning.

Recall Solution L5.1

Frequency: $ω_{s} = 9.8/6.371 \times 1 0^{6} = 1.2404 \times 1 0^{- 3} rad/s$ . Period: $T_{s} = 2 π / ω_{s} = 5065.6 s = 84.4 min$ — the classic Schuler period. Short-time limit: for $t ≪ T_{s}$ , the closed-loop solution behaves like $sin (ω_{s} t) \approx ω_{s} t - \frac{1}{6} (ω_{s} t)^{3}$ . Expanding the oscillatory closed-loop response to its lowest orders in $t$ reproduces exactly the polynomial (quadratic, then cubic) growth of the open-loop model. So the open-loop cubic law of L4.1 is simply the short-time limit of the bounded Schuler oscillation — it is the first bend of a swing that, given enough time, curves back rather than diverging. In one sentence: the $t^{3}$ growth is what the Schuler sine looks like before the spring has had a chance to pull the error home.

Recall Solution L5.2

Coupling: $- [f^{n} \times] ψ = 0$ when $f^{n} = 0$ , so tilt stops leaking into velocity during true free-fall. But not safe: the gyros keep drifting, so $\dot{ψ} = - C_{b}^{n} δ ω$ still grows the tilt $ψ$ . The moment specific force returns (engine burn, re-entry drag), that accumulated tilt suddenly re-projects a large force into velocity. Degenerate input silences the symptom, not the cause.

Recall Solution L5.3

Solve the inequality: $\frac{1}{6} g ϵ t^{3} \leq 50 \Rightarrow ϵ \leq \frac{6 \times 50}{9.8 \times 10 0 ^{3}} = \frac{300}{9.8 \times 1 0 ^{6}} = 3.0612 \times 1 0^{- 5} rad/s$ . Convert to $^{\circ} / hr$ : multiply by $\frac{180}{π}$ (rad→deg) and by $3600$ (per s → per hr): $ϵ \leq 3.0612 \times 1 0^{- 5} \times \frac{180}{π} \times 3600 = 6.314^{\circ} / hr$ . Design read: you need a gyro of $≲ 6^{\circ} / hr$ drift — comfortably a low-cost MEMS-to-tactical grade for a $100$ s flight. Longer flights ( $t^{3}$ !) tighten this brutally, which is why long-endurance INS demands navigation-grade gyros. This ties directly to Gyro and Accelerometer Error Models and the estimation limits in Kalman Filter for INS-GPS Integration.

Recall One-line summary to self-test

Count integrations from source to state ::: accel bias → position is quadratic ( $\frac{1}{2} b_{a} t^{2}$ ); gyro drift → position is cubic ( $\frac{1}{6} g ϵ t^{3}$ ); both are the short-time limbs of an $\approx 84$ min Schuler oscillation.