Estimating wave height from IMU accelerometer on a moving boat (Early Draft)

Usual methods for measuring elevation built-in into modern smartwatches and smart devices are based on fusion of data (using Kalman filters, etc) from IMU (accelerometer) and barometer sensors. These methods are designed to measure steps up and down stairs, hiking activities, and so on. However, using pressure as a useful source of elevation data is not suitable to measure the wave height from a boat. The surface of the wave itself (at least in trochoidal wave model) is the surface of equal pressure. Real waves are more complicated but nevertheless measuring air pressure on their surface to estimate amplitude of the wave would lead to measuring mostly noise.

Calculating elevation (displacement) using double integration of the acceleration leads to huge result drift due to not knowing initial displacement, initial vertical velocity and due to signal noise and accelerometer bias.

This article describes two different methods to estimate wave height, frequency, period from a moving boat using IMU (accelerometer) and other common boat sensors such as anemometer (wind), GPS, and compass. One method is based on Kalman method which calculates second integral of acceleration with a drift correction based on the fact that in wave average displacement equals to zero. Another method doesn’t use integration at all. Instead, it measures max/min acceleration and assumes wave shape to be trochoidal. Based on max/min acceleration and observed wave frequency (later adjusted using Doppler effect formulas) the method recreates parameters of the trochoid wave allowing it to estimate the wave height.

Initial velocity and position for Kalman filter are taken as estimates from trochoidal model to allow Kalman filter to converge faster.

Gravitational   \$g\$   (m/s^2):     \$g = 9.81\$

Rotation center of trochoidal wave particle in Z axis (always negative) (m):     \$b\$

Wavelength (m):     \$L\$

Wave equation:

\$X(t) = H sin(k c t)\$

\$Z(t) = - H cos(k c t)\$

Wavenumber (1/m):     \$k = (2 pi) / L\$

Wave amplitude (m):     \$H = e ^(k b) / k = L / (2 pi) e ^ ((2 pi b) / L) \$     ⇒     \$b = ln(H k) / k = L / (2 pi) ln ((2 pi H) / L) \$

max(H) amplitude (m):     \$H_max = L / (2 pi)\$     when   b = 0

Depth (m):     \$d\$

Wave speed (m/s):     \$c = sqrt(g / k tanh(k d))\$  ,      for deep water waves: \$c = sqrt(g / k) = sqrt((g L) / (2 pi))\$,    \$L = (g T ^ 2) / (2 pi)\$

Wave period (s):     \$T = L / c\$

Wave frequency (Hz=1/s):     \$f = 1 / T = c / L\$

Vertical motion is harmonic.

Vertical displacement: \$Z(t) = - H cos(k c t) = - L / (2 pi) (a(t))/g\$

Vertical velocity: \$u(t) = Z^'(t) = k c H sin(k c t) = c e ^ ((2 pi b) / L) sin(k c t)\$

\$u_max = c e ^ ((2 pi b) / L)\$

\$u_max(b=0) = c\$

Vertical acceleration: \$a(t) = u^'(t) = k^2 c^2 H cos(k c t) = g e ^ ((2 pi b) / L) cos(k c t)\$

\$a_max = -a_min = k^2 c^2 H = ((2 pi) / L) ^ 2 ((g L) / (2 pi)) L / (2 pi) e^((2 pi b)/L) = g e^((2 pi b)/L) \$

\$b = L / (2 pi) ln(a_max / g)\$

Slope \$s(t)\$:

\$x_b(t) = X(t) - t Delta v\$

\$s(t) = d/ (dx_b) ((Z(t)) / (x_b(t))) = (d/dt(Z)) / (d/dt(x_b)) = (d/dt(- H cos(kct))) / (d/dt(H sin(kct) - t Delta v)) = (cHk sin(kct)) / (cHk cos(kct) - Delta v)\$

Ref:

Sharkh S. M., Hendijanizadeh2 M., Moshrefi-Torbati3 M., Abusara M. A.: A Novel Kalman Filter Based Technique for Calculating the Time History of Vertical Displacement of a Boat from Measured Acceleration, Marine Engineering Frontiers (MEF) Volume 2, 2014

The method uses assumption that average vertical displacement is zero. So it takes third integral of acceleration as input measurement with zero value. Vertical acceleration on each step is the value used in transition offset.

So only input observation value is 0 (third integral of acceleration). Acceleration used in transition offset needs to exclude Earth gravitational g and accelerometer bias. So avg(acceleration) is subtracted from it.

Accelerometer noise will make the filter to converge slower. Accelerometer bias will skew the displacement on one side. Giving too low covariances will flatten the result too much. Too high covariances will make it jump too much.

Low pass filter is needed to get rid of high frequency noise before running this filter to get better results.

Kalman method parameters.

Transition matrix:

\$ F = [[1, dt, 1/2 dt ^ 2 ], [0, 1, dt], [0, 0, 1]] \$

Transition offset:

\$ B = ( (1/6 dt^3), (1/2 dt^2), (dt) ) \$

Observation matrix:

\$ H = (1, 0, 0) \$

Transition covariance:

\$ Q = [[ \text{PosIntegral_Trans_Variance}, 0, 0 ], [0, \text{Pos_Trans_Variance}, 0 ], [0, 0, \text{Vel_Trans_Variance} ]] \$

Observation covariance:

\$ R = [[ \text{PosIntegral_Variance} ]] \$

Initial state mean (height integral, height, velocity):

\$ X0 = ((0), (\text{height_from_trochoid_model}), (\text{velocity_from_trochoid_model})) \$

Initial state covariance:

\$ P0 = [[ \text{PosIntegral_Variance}, 0, 0], [0, 1, 0], [0, 0, 1]] \$

The wave height needs to be evaluated based on the results on the end side of the sampling interval. The results in the beginning of the sample interval will have extra spikes due to Kalman filter taking time to converge (without knowing initial velocity and position and adjusting covariances).

Trochoidal model gives a good reference point. It’s easy to go from acceleration to vertical displacement in trochoidal model knowing the wavelength. So (I think) it could be a good idea to use trochoidal model estimations as input to Kalman filter fusion (instead of zero third integral assumption), as another way to prevent displacement integral of acceleration from drift.

Vertical velocity and heel are correlated. It would be interesting to study their interdependencies.