Marine IMU Wave-State Estimation with OU-Driven Quaternion MEKF

This document replaces the earlier trochoidal-wave and raw double-integration draft with the more mature estimator now developed in the Ocean IMU project. The current method is an OU-driven quaternion multiplicative extended Kalman filter, or Q-MEKF, for marine inertial sensing and wave-state estimation. It estimates attitude, heave, 3D displacement, dominant wave frequency, and horizontal wave direction from a calibrated IMU, with optional magnetometer yaw correction and vessel-specific corrections such as IMU lever arm and steady wind heel.

The central change is conceptual. The old draft tried to infer wave height from simplified trochoidal formulas and then use those formulas to help a vertical Kalman integration scheme. The current method instead models the unknown world-frame wave acceleration as a stationary Ornstein-Uhlenbeck, or OU, process and lets the Kalman filter jointly estimate attitude and wave kinematics. Drift is controlled by a soft pseudo-measurement on the integral of displacement, with an adaptive covariance that scales with the sea state.

This formulation is intended for embedded, low-cost marine IMU deployments. It avoids reliance on barometric pressure for wave elevation, avoids uncontrolled double integration of acceleration, and uses a real-time adaptation loop so that the same filter can remain stable across calm chop, long swell, and larger sea states.

This document is a practical AsciiDoc distillation of the more detailed Ocean IMU mathematical draft:

It is not meant to duplicate every proof, coefficient table, or plotting fixture from Ocean IMU. Instead, it records the model, state layout, update equations, tuning law, and validation results at a level appropriate for the Bareboat Necessities documentation.

A boat, buoy, or floating sensor experiences wave-induced translational motion, vessel attitude motion, sensor bias, vibration, magnetic disturbance, and measurement noise at the same time. Estimating wave height by simply rotating the accelerometer into a vertical axis and double-integrating it is not stable in practice:

The mature estimator treats the problem as a coupled inertial-navigation and stochastic wave-kinematics problem. The wave acceleration is not assumed to be a perfect sinusoid or trochoid. Instead, it is treated as colored, finite-variance excitation with a correlation time tied to the dominant wave period.

The estimator combines four pieces:

At each IMU step:

The working world frame is NED: north, east, down. The body frame is the calibrated IMU or vessel body frame used by the implementation.

The estimator assumes calibrated triads:

The gravity vector in NED is

The rotation \$R_wb\$ maps world-frame vectors into the body frame. The nominal attitude is stored as a unit quaternion representing world-to-body orientation.

The extended error-state vector is

The state components are:

With all components enabled, the filter has

states.

For sample interval \$h\$, a linear time-invariant SDE is discretized as

The Ocean IMU implementation computes the important OU-driven blocks analytically. This avoids numerical instability and avoids repeatedly using a general-purpose matrix exponential on small embedded targets.

The accelerometer, magnetometer, and pseudo-measurement channels are all rank-3 updates. For a measurement

the EKF update is

The covariance update should use the Joseph form:

After each correction, the attitude-error subvector is injected back into the nominal quaternion:

The rank-3 structure is important for embedded use. The expensive inversion is only a 3 by 3 solve for each triad update.

Startup should be staged. The filter should not immediately trust yaw, acceleration bias, displacement, or drift-control states.

A practical sequence is:

This avoids seeding displacement, bias, and pseudo-measurement states with startup transients.

Adaptive tuning requires an online estimate of the dominant wave frequency. The signal is derived from vertical inertial acceleration, typically as an up-positive scalar:

Three embedded-friendly trackers are supported in Ocean IMU:

The raw frequency is clamped and smoothed to produce \$f_tune\$, which drives the sea-state tuner and the horizontal wave-direction estimator.

The tuner adapts three quantities:

The variance estimate is tied to a fixed number of wave periods so that the adaptation horizon scales naturally between chop and swell:

After subtracting a configured acceleration noise floor,

the target parameters are

The applied values are exponentially smoothed. \$tau\$ and \$sigma_a\$ use a fixed adaptation time constant. \$R_S\$ is smoothed with a time constant proportional to the current wave time scale, so it responds faster in short-period seas and slower in long-period swell.

The filter estimates a slowly varying horizontal acceleration amplitude vector \$A_k\$. Let

For a dominant narrow-band wave with phase \$phi_k\$, the simplified observation model is

This is the online Kalman-filter form of a cosine-projection least-squares fit:

The 2D direction filter uses

When \$abs(cos(phi_k))\$ is too small, the update is skipped because the observation carries little information near horizontal zero crossings.

The wave line angle is

This gives direction modulo \$180 deg\$. A separate vertical-horizontal phase discriminator can choose the travel sense along that line, producing a directed azimuth when the sign estimate is stable.

Ocean IMU validates the filter on synthetic directional sea states with controlled ground truth. The simulations use 20 minute runs at 200 Hz, with IMU and magnetometer error models representative of low-cost MEMS hardware.

The implementation has also been exercised on real embedded hardware using the M5Stack AtomS3R ESP32-S3 module. The purpose of this validation is to confirm real-time behavior and implementation feasibility:

Quantitative wave-kinematics accuracy is evaluated in simulation, where exact displacement, velocity, acceleration, and attitude ground truth are available. Real hardware validates the embedded execution path and sensor-processing pipeline.

The old document contained useful early notes about why double integration drifts and why pressure is not a reliable wave-height source on the wave surface. Those points remain valid. However, the trochoidal reconstruction and Doppler-based approach is no longer the preferred core algorithm.

The current preferred algorithm is: