Master Thesis Code
by Simon Moser
Loading...
Searching...
No Matches
Part 7: Initial Parameter Estimation

Last updated: 22. April 2024

This section introduces the considerations of an optimization algorithm for initial parameter estimation. First, the measurement model is defined. Then, a nonlinear least-suares optimization algorithm is used to estimate the calibration parameters of the system.

Measurement Model

The measurement models are based on Guang-Lin et al., (2012).

Acceleration

The acceleration measurement at a constant temperature can be modeled as

\begin{align*} a_x^\text{(meas)} &= b_x^{(a)} + s_{xx}^{(a)} a_{x} + s_{xy}^{(a)} a_{y} + s_{xz}^{(a)} a_{z} + s_{\omega_{xx}}^{(a)} \omega_{x} + s_{\omega_{xy}}^{(a)} \omega_{y} + s_{\omega_{xz}^{(a)}} \omega_{z} + \varepsilon_{ax} \\[3mm] a_y^\text{(meas)} &= b_y^{(a)} + s_{yx}^{(a)} a_{x} + s_{yy}^{(a)} a_{y} + s_{yz}^{(a)} a_{z} + s_{\omega_{yx}}^{(a)} \omega_{x} + s_{\omega_{yy}}^{(a)} \omega_{y} + s_{\omega_{yz}^{(a)}} \omega_{z} + \varepsilon_{ay} \\[3mm] a_z^\text{(meas)} &= b_z^{(a)} + s_{zx}^{(a)} a_{x} + s_{zy}^{(a)} a_{y} + s_{zz}^{(a)} a_{z} + s_{\omega_{zx}}^{(a)} \omega_{x} + s_{\omega_{zy}}^{(a)} \omega_{y} + s_{\omega_{zz}^{(a)}} \omega_{z} + \varepsilon_{az}, \end{align*}

where

  • \(a_i^\text{(meas)}\) is the measured acceleration in the \(i\)-th direction,
  • \(b_i^\text{(a)}\) is the constant bias of the accelerometer in the \(i\)-th direction,
  • \(s_{ii}^{(a)}\) is the scale factor of the accelerometer in the \(i\)-th direction,
  • \(s_{ij}^{(a)}\) is the cross-axis sensitivity of the accelerometer in the \(i\)-th direction to the \(j\)-th direction,
  • \(s_{\omega_{ij}}^{(a)}\) is the cross-sensor sensitivity of the gyroscope in the \(i\)-th direction to the \(j\)-th direction,
  • \(\varepsilon_{ai}\) is the residual error of the accelerometer in the \(i\)-th direction (e.g. due to temperature changes or random processes),

with \(i,j \in \{x,y,z\}\).

Gyroscope

The gyroscope measurement at a constant temperature can be modeled as

\begin{align*} \omega_x^\text{(meas)} &= b_x^{(\omega)} + s_{xx}^{(\omega)} \omega_{x} + s_{xy}^{(\omega)} \omega_{y} + s_{xz}^{(\omega)} \omega_{z} + s_{a_{xx}}^{(\omega)} a_{x} + s_{a_{xy}}^{(\omega)} a_{y} + s_{a_{xz}}^{(\omega)} a_{z} + s_{a_{2x}}^{(\omega)} a_x^2 + \varepsilon_{\omega x} \\[3mm] \omega_y^\text{(meas)} &= b_y^{(\omega)} + s_{yx}^{(\omega)} \omega_{x} + s_{yy}^{(\omega)} \omega_{y} + s_{yz}^{(\omega)} \omega_{z} + s_{a_{yx}}^{(\omega)} a_{x} + s_{a_{yy}}^{(\omega)} a_{y} + s_{a_{yz}}^{(\omega)} a_{z} + s_{a_{2y}}^{(\omega)} a_y^2 + \varepsilon_{\omega y} \\[3mm] \omega_z^\text{(meas)} &= b_z^{(\omega)} + s_{zx}^{(\omega)} \omega_{x} + s_{zy}^{(\omega)} \omega_{y} + s_{zz}^{(\omega)} \omega_{z} + s_{a_{zx}}^{(\omega)} a_{x} + s_{a_{zy}}^{(\omega)} a_{y} + s_{a_{zz}}^{(\omega)} a_{z} + s_{a_{2z}}^{(\omega)} a_z^2 + \varepsilon_{\omega z}, \end{align*}

where

  • \(\omega_i^\text{(meas)}\) is the measured angular velocity in the \(i\)-th direction,
  • \(b_i^{(\omega)}\) is the constant bias of the gyroscope in the \(i\)-th direction,
  • \(s_{ii}^{(\omega)}\) is the scale factor of the gyroscope in the \(i\)-th direction,
  • \(s_{ij}^{(\omega)}\) is the cross-axis sensitivity of the gyroscope in the \(i\)-th direction to the \(j\)-th direction,
  • \(s_{a_{ij}}^{(\omega)}\) is the cross-sensor sensitivity of the accelerometer in the \( i \)-th direction to the \( j \)-th direction,
  • \(s_{a_{2i}}^{(\omega)}\) is the cross-sensor sensitivity of the accelerometer in the \( i \)-th direction to the \( j \)-th direction,
  • \(\varepsilon_{\omega i}\) is the residual error of the gyroscope in the \(i\)-th direction (e.g. due to temperature changes or random processes),

with \(i,j \in \{x,y,z\}\).

Simplified Measurement Model 1

If we only consider the accelerometer and gyroscope measurements during zero velocity periods, we know that

\begin{align*} \omega_i &= 0, \quad \forall \ i \in \{x,y,z\} \\ || a ||_2 &= g, \quad \text{where } g \text{ is the gravitational acceleration}. \end{align*}

Note
With the zero velocity assumption, the sensitivity of the gyroscope cannot be determined.

Thus, the measurement models can be simplified to

\begin{align*} a_x^\text{(meas)} &= b_x^{(a)} + s_{xx}^{(a)} a_{x} + s_{xy}^{(a)} a_{y} + s_{xz}^{(a)} a_{z} + \varepsilon_{ax} \\[3mm] a_y^\text{(meas)} &= b_y^{(a)} + s_{yx}^{(a)} a_{x} + s_{yy}^{(a)} a_{y} + s_{yz}^{(a)} a_{z} + \varepsilon_{ay} \\[3mm] a_z^\text{(meas)} &= b_z^{(a)} + s_{zx}^{(a)} a_{x} + s_{zy}^{(a)} a_{y} + s_{zz}^{(a)} a_{z} + \varepsilon_{az}, \end{align*}

and

\begin{align*} \omega_x^\text{(meas)} &= b_x^{(\omega)} + s_{a_{xx}}^{(\omega)} a_{x} + s_{a_{xy}}^{(\omega)} a_{y} + s_{a_{xz}}^{(\omega)} a_{z} + s_{a_{2x}}^{(\omega)} a_x^2 + \varepsilon_{\omega x} \\[3mm] \omega_y^\text{(meas)} &= b_y^{(\omega)} + s_{a_{yx}}^{(\omega)} a_{x} + s_{a_{yy}}^{(\omega)} a_{y} + s_{a_{yz}}^{(\omega)} a_{z} + s_{a_{2y}}^{(\omega)} a_y^2 + \varepsilon_{\omega y} \\[3mm] \omega_z^\text{(meas)} &= b_z^{(\omega)} + s_{a_{zx}}^{(\omega)} a_{x} + s_{a_{zy}}^{(\omega)} a_{y} + s_{a_{zz}}^{(\omega)} a_{z} + s_{a_{2z}}^{(\omega)} a_z^2 + \varepsilon_{\omega z}. \end{align*}

Even tough these equations are simplified, there are still twelve or fifeteen unknown parameters, respevtively. Therefore, a calibration sequence of at least fifeen different known poses is required to estimate all parameters.

Simplified Measurement Model 2

A duable calibration sequence is to put the sensor on a flat surface where every face is pointing upwards once. If we consider such a calibration sequence where the sensor is in six different poses, where always one axis is aligned with the gravitational acceleration, the calibration model must be further simplified.

Since for our sensors (MPU-9250 and ADXL375) the cross-axis sensitivities are assumed to be smaller than the scale factors, we neglect the cross-axis sensitivities. Furthermore, we assume that the gyroscope is not very sensitive to the accelerometer measurements (small cross-sensor sensitivities). Thus, the measurement models can be further simplified to

\begin{align*} a_x^\text{(meas)} &= b_x^{(a)} + s_{xx}^{(a)} a_{x} + \varepsilon_{ax} \\[3mm] a_y^\text{(meas)} &= b_y^{(a)} + s_{yy}^{(a)} a_{y} + \varepsilon_{ay} \\[3mm] a_z^\text{(meas)} &= b_z^{(a)} + s_{zz}^{(a)} a_{z} + \varepsilon_{az} \end{align*}

for the accelerometer and

\begin{align*} \omega_x^\text{(meas)} &= b_x^{(\omega)} + \varepsilon_{\omega x} \\[3mm] \omega_y^\text{(meas)} &= b_y^{(\omega)} + \varepsilon_{\omega y} \\[3mm] \omega_z^\text{(meas)} &= b_z^{(\omega)} + \varepsilon_{\omega z} \end{align*}

for the gyroscope.

Algorithm to estimate the calibration parameters

A function jumpEstimateCalibrationParameters.m is implemented to estimate the calibration parameters. This function takes sensor data and a zero velocity detection vector as inputs. The general functionality is as follows:

  • Detect the start and end indices of the zero velocity periods based on the zero velocity detection vector.
  • Remove zero velocity periods that are shorter than a threshold (option minZvDuration).
  • Split the zero velocity periods into smaller segments if the length of the zero velocity period is larger than a threshold (option maxZvDuration).
  • Calculate the mean of the accelerometer and high-g accelerometer measurements for each segment.
  • Detect the orientation of the sensor based on the mean accelerometer measurements, since we know that one of the axes is aligned with the gravitational acceleration.
  • Detect zero-velocity periods where the sensor is not aligned with the gravitational acceleration and remove them.
  • Check that every of the six possible orientations is detected at least once.
  • Estimate the calibration parameters based on the detected zero velocity periods by minimizing the residuals \(\varepsilon_{ij}\) for each sensor measurement.

The function jumpEstimateCalibrationParameters.m is based on the MATLAB function lsqlin to solve the optimization problem.

Example Implementation

An example of the usage of the function is shown in derivation7IntialParamEst.m. There, the parameters are estimated three times:

  1. Using synthetic data with a known ground truth.
  2. Using actual sensor data from a calibration sequence.
  3. Using actual sensor data from a calibration sequence where some of the zero velocity periods were intentionally not aligned with the gravitational acceleration (to test the robustness of the algorithm).

The output of the text in the example is shown below.

Sequence 1, (synthetic data)
MPU Accelerometer Parameters:
Known Parameters from simulation:
Biases: x=-0.382, y=0.695, z=-0.595
Sensitivities: x=0.967, y=0.975, z=0.983
Estimated Parameters from data:
Biases: x=-0.383, y=0.694, z=-0.594
Sensitivities: x=0.967, y=0.974, z=0.983
High G Accelerometer Parameters:
Known Parameters from simulation:
Biases: x=2.368, y=-2.100, z=-0.609
Sensitivities: x=1.000, y=1.000, z=1.000
Estimated Parameters from data:
Biases: x=2.373, y=-2.096, z=-0.628
Sensitivities: x=1.000, y=1.001, z=1.002
=======================================
Sequence 2, (real data)
Estimated Accel Parameters from real data:
Biases: x=0.255, y=0.115, z=0.760
Sensitivities: x=1.000, y=1.001, z=0.999
Estimated High G Accel Parameters from real data:
Biases: x=-7.946, y=0.185, z=-1.146
Sensitivities: x=0.942, y=0.970, z=0.995
=======================================
Warning: Found invalid range and removed it. Was the sensor laying flat on a surface?
> In jumpEstimateCalibrationParameters (line 98)
In derivation7IntialParamEst (line 83)
Warning: Found invalid range and removed it. Was the sensor laying flat on a surface?
> In jumpEstimateCalibrationParameters (line 98)
In derivation7IntialParamEst (line 83)
Warning: Found invalid range and removed it. Was the sensor laying flat on a surface?
> In jumpEstimateCalibrationParameters (line 98)
In derivation7IntialParamEst (line 83)
Sequence 3, (real data, not always aligned)
Estimated Parameters from real data with non-flat cube:
Biases: x=0.271, y=0.135, z=0.784
Sensitivities: x=1.000, y=1.001, z=0.999
Estimated High G Accel Parameters from real data with non-flat cube:
Biases: x=-7.562, y=0.273, z=-1.341
Sensitivities: x=0.946, y=0.965, z=0.987