insfilterErrorState

Estimate pose from IMU, GPS, and monocular visual odometry (MVO) data

Description

The insfilterErrorState object implements sensor fusion of IMU, GPS, and monocular visual odometry (MVO) data to estimate pose in the NED (or ENU) reference frame. The filter uses a 17-element state vector to track the orientation quaternion, velocity, position, IMU sensor biases, and the MVO scaling factor. The insfilterErrorState object uses an error-state Kalman filter to estimate these quantities.

Creation

Syntax

filter = insfilterErrorState

filter = insfilterErrorState('ReferenceFrame',RF)

filter = insfilterErrorState(___,Name=Value)

Description

filter = insfilterErrorState creates an insfilterErrorState object with default property values.

filter = insfilterErrorState('ReferenceFrame',RF) allows you to specify the reference frame, RF, of the filter.

filter = insfilterErrorState(___,Name=Value) sets one or more properties using name-value arguments in addition to any of the previous input arguments.

Input Arguments

RF — Local navigation coordinate system

'NED' (default) | 'ENU'

Local navigation coordinate system, specified as either 'NED' (North-East-Down) or 'ENU' (East-North-Up).

Data Types: char | string

Properties

IMUSampleRate — Sample rate of IMU (Hz)

100 (default) | positive scalar

Sample rate of the inertial measurement unit (IMU) in Hz, specified as a positive scalar.

Data Types: single | double

ReferenceLocation — Reference location (deg, deg, meters)

[0 0 0] (default) | 3-element positive row vector

Reference location, specified as a 3-element row vector in geodetic coordinates (latitude, longitude, and altitude). Altitude is the height above the reference ellipsoid model, WGS84. The reference location units are [degrees degrees meters].

Data Types: single | double

GyroscopeNoise — Multiplicative process noise variance from gyroscope ((rad/s)²)

[1e-6 1e-6 1e-6] (default) | scalar | 3-element row vector

Multiplicative process noise variance from the gyroscope in (rad/s)², specified as a scalar or 3-element row vector of positive real finite numbers.

• If GyroscopeNoise is specified as a row vector, the elements correspond to the noise in the x, y, and z axes of the gyroscope, respectively.

• If GyroscopeNoise is specified as a scalar, the single element is applied to each axis.

Data Types: single | double

GyroscopeBiasNoise — Additive process noise variance from gyroscope bias ((rad/s)²)

[1e-9 1e-9 1e-9] (default) | scalar | 3-element row vector

Additive process noise variance from the gyroscope bias in (rad/s)², specified as a scalar or 3-element row vector of positive real finite numbers.

• If GyroscopeBiasNoise is specified as a row vector, the elements correspond to the noise in the x, y, and z axes of the gyroscope, respectively.

• If GyroscopeBiasNoise is specified as a scalar, the single element is applied to each axis

Data Types: single | double

AccelerometerNoise — Multiplicative process noise variance from accelerometer ((m/s²)²)

[1e-4 1e-4 1e-4] (default) | scalar | 3-element row vector

Multiplicative process noise variance from the accelerometer in (m/s²)², specified as a scalar or 3-element row vector of positive real finite numbers.

• If AccelerometerNoise is specified as a row vector, the elements correspond to the noise in the x, y, and z axes of the accelerometer, respectively.

• If AccelerometerNoise is specified as a scalar, the single element is applied to each axis.

Data Types: single | double

AccelerometerBiasNoise — Additive process noise variance from accelerometer bias ((m/s²)²)

[1e-4 1e-4 1e-4] (default) | scalar | 3-element row vector

Additive process noise variance from accelerometer bias in (m/s²)², specified as a scalar or 3-element row vector of positive real numbers.

• If AccelerometerBiasNoise is specified as a row vector, the elements correspond to the noise in the x, y, and z axes of the accelerometer, respectively.

• If AccelerometerBiasNoise is specified as a scalar, the single element is applied to each axis.

State — State vector of Kalman filter

[1;zeros(15,1);1] (default) | 17-element column vector

State vector of the extended Kalman filter, specified as a 17-element column vector. The state values represent:

State	Units	Index
Orientation (quaternion parts)	N/A	1:4
Position (NED or ENU)	m	5:7
Velocity (NED or ENU)	m/s	8:10
Gyroscope Bias (XYZ)	rad/s	11:13
Accelerometer Bias (XYZ)	m/s2	14:16
Visual Odometry Scale (XYZ)	N/A	17

The default initial state corresponds to an object at rest located at [0 0 0] in geodetic LLA coordinates.

Data Types: single | double

StateCovariance — State error covariance for Kalman filter

ones(16) (default) | 16-by-16 matrix

State error covariance for the Kalman filter, specified as a 16-by-16-element matrix of real numbers. The state error covariance values represent:

State Covariance	Row/Column Index
δ Rotation Vector (XYZ)	1:3
δ Position (NED or ENU)	4:6
δ Velocity (NED or ENU)	7:9
δ Gyroscope Bias (XYZ)	10:12
δ Accelerometer Bias (XYZ)	13:15
δ Visual Odometry Scale (XYZ)	16

Note that because this is an error-state Kalman filter, it tracks the errors in the states. δ represents the error in the corresponding state.

Data Types: single | double

Object Functions

predict	Update states using accelerometer and gyroscope data for insfilterErrorState
correct	Correct states using direct state measurements for insfilterErrorState
residual	Residuals and residual covariances from direct state measurements for insfilterErrorState
fusegps	Correct states using GPS data for insfilterErrorState
residualgps	Residuals and residual covariance from GPS measurements for insfilterErrorState
fusemvo	Correct states using monocular visual odometry for insfilterErrorState
residualmvo	Residuals and residual covariance from monocular visual odometry measurements for insfilterErrorState
pose	Current orientation and position estimate for insfilterErrorState
reset	Reset internal states for insfilterErrorState
stateinfo	Display state vector information for insfilterErrorState
tune	Tune insfilterErrorState parameters to reduce estimation error
copy	Create copy of insfilterErrorState

Examples

Estimate Pose of Ground Vehicle

Load logged data of a ground vehicle following a circular trajectory. The .mat file contains IMU and GPS sensor measurements and ground truth orientation and position.

load('loggedGroundVehicleCircle.mat', ...
    'imuFs','localOrigin', ...
    'initialStateCovariance', ...
    'accelData','gyroData', ...
    'gpsFs','gpsLLA','Rpos','gpsVel','Rvel', ...
    'trueOrient','truePos');

Create an INS filter to fuse IMU and GPS data using an error-state Kalman filter.

initialState = [compact(trueOrient(1)),truePos(1,:),-6.8e-3,2.5002,0,zeros(1,6),1].';
filt = insfilterErrorState;
filt.IMUSampleRate = imuFs;
filt.ReferenceLocation = localOrigin;
filt.State = initialState;
filt.StateCovariance = initialStateCovariance;

Preallocate variables for position and orientation. Allocate a variable for indexing into the GPS data.

numIMUSamples = size(accelData,1);
estOrient = ones(numIMUSamples,1,'quaternion');
estPos = zeros(numIMUSamples,3);

gpsIdx = 1;

Fuse accelerometer, gyroscope, and GPS data. The outer loop predicts the filter forward at the fastest sample rate (the IMU sample rate).

for idx = 1:numIMUSamples

    % Use predict to estimate the filter state based on the accelData and
    % gyroData arrays.
    predict(filt,accelData(idx,:),gyroData(idx,:));
    
    % GPS data is collected at a lower sample rate than IMU data. Fuse GPS
    % data at the lower rate.
    if mod(idx, imuFs / gpsFs) == 0
        % Correct the filter states based on the GPS data.
        fusegps(filt,gpsLLA(gpsIdx,:),Rpos,gpsVel(gpsIdx,:),Rvel);
        gpsIdx = gpsIdx + 1;
    end
    
    % Log the current pose estimate
    [estPos(idx,:), estOrient(idx,:)] = pose(filt);
end

Calculate the RMS errors between the known true position and orientation and the output from the error-state filter.

pErr = truePos - estPos;
qErr = rad2deg(dist(estOrient,trueOrient));

pRMS = sqrt(mean(pErr.^2));
qRMS = sqrt(mean(qErr.^2));

fprintf('Position RMS Error\n');

Position RMS Error

fprintf('\tX: %.2f, Y: %.2f, Z: %.2f (meters)\n\n',pRMS(1),pRMS(2),pRMS(3));

	X: 0.40, Y: 0.24, Z: 0.05 (meters)

fprintf('Quaternion Distance RMS Error\n');

Quaternion Distance RMS Error

fprintf('\t%.2f (degrees)\n\n',qRMS);

	0.30 (degrees)

Visualize the true position and the estimated position.

plot(truePos(:,1),truePos(:,2),estPos(:,1),estPos(:,2),'r:','LineWidth',2)
grid on
axis square
xlabel('N (m)')
ylabel('E (m)')
legend('Ground Truth','Estimation')

Algorithms

Note: The following algorithm only applies to an NED reference frame.

insfilterErrorState uses a 17-axis error state Kalman filter structure to estimate pose in the NED reference frame. The state is defined as:

$$x=\left[\begin{array}{c}{q}_0\\ q_1\\ q_2\\ q_3\\ positio{n}_N\\ positio{n}_E\\ positio{n}_D\\ v_N\\ v_E\\ v_D\\ gyrobia{s}_X\\ gyrobia{s}_Y\\ gyrobia{s}_Z\\ accelbia{s}_X\\ accelbia{s}_Y\\ \begin{array}{l}accelbia{s}_Z\\ \begin{array}{c}scaleFactor\end{array}\end{array}\end{array}\right]$$

where

• q₀, q₁, q₂, q₃ –– Parts of orientation quaternion. The orientation quaternion represents a frame rotation from the platform's current orientation to the local NED coordinate system.

• position_N, position_E, position_D –– Position of the platform in the local NED coordinate system.

• gyrobias_X, gyrobias_Y, gyrobias_Z –– Bias in the gyroscope reading.

• accelbias_X, accelbias_Y, accelbias_Z –– Bias in the accelerometer reading.

• scaleFactor –– Scale factor of the pose estimate.

Given the conventional formulation of the state transition function,

$$x_{k|k-1}=f({\widehat{x}}_{k-1|k-1})$$

the predicted state estimate is:

$$x_{k|k-1}=\left[\begin{array}{c}{q}_0+\Delta t\ast q_1(gyrobia{s}_X/2-gyr{o}_X/2)+\Delta t\ast q_2\ast (gyrobia{s}_Y/2-gyr{o}_Y/2)+\Delta t\ast q_3\ast (gyrobia{s}_Z/2-gyr{o}_Z/2)\\ q_1-\Delta t\ast q_0(gyrobia{s}_X/2-gyr{o}_X/2)+\Delta t\ast q_3\ast (gyrobia{s}_Y/2-gyr{o}_Y/2)-\Delta t\ast q_2\ast (gyrobia{s}_Z/2-gyr{o}_Z/2)\\ q_2-\Delta t\ast q_3(gyrobia{s}_X/2-gyr{o}_X/2)-\Delta t\ast q_0\ast (gyrobia{s}_Y/2-gyr{o}_Y/2)+\Delta t\ast q_1\ast (gyrobia{s}_Z/2-gyr{o}_Z/2)\\ q_3+\Delta t\ast q_2(gyrobia{s}_X/2-gyr{o}_X/2)-\Delta t\ast q_1\ast (gyrobia{s}_Y/2-gyr{o}_Y/2)-\Delta t\ast q_0\ast (gyrobia{s}_Z/2-gyr{o}_Z/2)\\ positio{n}_N+\Delta t\ast v_N\\ positio{n}_E+\Delta t\ast v_E\\ positio{n}_D+\Delta t\ast v_D\\ v_N-\Delta t\ast \left\{\begin{array}{l}{q}_0\ast \left(q_0\ast \left(accelbia{s}_X-acce{l}_X\right)-q_3\ast \left(accelbia{s}_Y-acce{l}_Y\right)+q_2\ast \left(accelbia{s}_Z-acce{l}_Z\right)\right)+g_N+\\ q_2\ast \left(q_1\ast \left(accelbia{s}_Y-acce{l}_Y\right)-q_2\ast \left(accelbia{s}_X-acce{l}_X\right)+q_0\ast \left(accelbia{s}_Z-acce{l}_Z\right)\right)+\\ q_1\ast \left(q_1\ast \left(accelbia{s}_X-acce{l}_X\right)+q_2\ast \left(accelbia{s}_Y-acce{l}_Y\right)+q_3\ast \left(accelbia{s}_Z-acce{l}_Z\right)\right)-\\ q_3\ast \left(q_3\ast \left(accelbia{s}_X-acce{l}_X\right)+q_0\ast \left(accelbia{s}_Y-acce{l}_Y\right)-q_1\ast \left(accelbia{s}_Z-acce{l}_Z\right)\right)\end{array}\right\}\\ v_E-\Delta t\ast \left\{\begin{array}{l}{q}_0\ast \left(q_3\ast \left(accelbia{s}_X-acce{l}_X\right)+q_0\ast \left(accelbia{s}_Y-acce{l}_Y\right)-q_1\ast \left(accelbia{s}_Z-acce{l}_Z\right)\right)+g_E-\\ q_1\ast \left(q_1\ast \left(accelbia{s}_Y-acce{l}_Y\right)-q_2\ast \left(accelbia{s}_X-acce{l}_X\right)+q_0\ast \left(accelbia{s}_Z-acce{l}_Z\right)\right)+\\ q_2\ast \left(q_1\ast \left(accelbia{s}_X-acce{l}_X\right)+q_2\ast \left(accelbia{s}_Y-acce{l}_Y\right)+q_3\ast \left(accelbia{s}_Z-acce{l}_Z\right)\right)+\\ q_3\ast \left(q_0\ast \left(accelbia{s}_X-acce{l}_X\right)-q_3\ast \left(accelbia{s}_Y-acce{l}_Y\right)+q_2\ast \left(accelbia{s}_Z-acce{l}_Z\right)\right)\end{array}\right\}\\ v_D-\Delta t\ast \left\{\begin{array}{l}{q}_0\ast \left(q_1\ast \left(accelbia{s}_Y-acce{l}_Y\right)-q_2\ast \left(accelbia{s}_X-acce{l}_X\right)+q_0\ast \left(accelbia{s}_Z-acce{l}_Z\right)\right)+g_D+\\ q_1\ast \left(q_3\ast \left(accelbia{s}_X-acce{l}_X\right)+q_0\ast \left(accelbia{s}_Y-acce{l}_Y\right)-q_1\ast \left(accelbia{s}_Z-acce{l}_Z\right)\right)-\\ q_2\ast \left(q_0\ast \left(accelbia{s}_X-acce{l}_X\right)-q_3\ast \left(accelbia{s}_Y-acce{l}_Y\right)+q_2\ast \left(accelbia{s}_Z-acce{l}_Z\right)\right)-\\ q_3\ast \left(q_1\ast \left(accelbia{s}_X-acce{l}_X\right)+q_2\ast \left(accelbia{s}_Y-acce{l}_Y\right)+q_3\ast \left(accelbia{s}_Z-acce{l}_Z\right)\right)\end{array}\right\}\\ gyrobia{s}_X\\ gyrobia{s}_Y\\ gyrobia{s}_Z\\ accelbia{s}_X\\ accelbia{s}_Y\\ accelbia{s}_Z\\ scaleFactor\end{array}\right]$$

where

• Δt –– IMU sample time.

• g_N, g_E, g_D –– Constant gravity vector in the NED frame.

Extended Capabilities

C/C++ Code Generation
Generate C and C++ code using MATLAB® Coder™.

Version History

Introduced in R2019a