ecompass

Ориентация от показаний магнитометра и акселерометра

Синтаксис

orientation = ecompass(accelerometerReading,magnetometerReading)
orientation = ecompass(accelerometerReading,magnetometerReading,orientationFormat)

Описание

пример

orientation = ecompass(accelerometerReading,magnetometerReading) возвращает quaternion, который может вращать количества от родительского элемента (NED) кадр дочернему элементу (датчик) кадр.

пример

orientation = ecompass(accelerometerReading,magnetometerReading,orientationFormat) задает формат ориентации как матрицу вращения или quaternion.

Примеры

свернуть все

Используйте известную силу магнитного поля, и соответствующее ускорение устройства указало истинный север в Бостоне, чтобы определить магнитный наклон Бостона.

Задайте известное ускорение и силу магнитного поля в Бостоне.

magneticFieldStrength = [19.535 -5.109 47.930];
properAcceleration = [0 0 9.8];

Передайте силу магнитного поля и ускорение к функции ecompass. Функция ecompass возвращает оператор вращения кватерниона. Преобразуйте кватернион в Углы Эйлера в градусах.

q = ecompass(properAcceleration,magneticFieldStrength);
e = eulerd(q,'ZYX','frame');

Угол, e, представляет угол между истинным северным и магнитным севером в Бостоне. Условно, магнитный наклон отрицателен, когда магнитный север к западу от истинного севера. Инвертируйте угол, чтобы определить магнитный наклон.

magneticDeclinationOfBoston = -e(1)
magneticDeclinationOfBoston =

  -14.6563

Функция ecompass плавит данные о магнитометре и акселерометре, чтобы возвратить кватернион, который, когда используется в операторе вращения кватерниона, может вращать количества от родительского элемента (NED) кадр к дочернему кадру. Функция ecompass может также возвратить матрицы вращения, которые выполняют эквивалентные вращения как оператор кватерниона.

Задайте вращение, которое может взять родительский кадр, указывающий на магнитный север к дочернему кадру, указывающему на географический север. Задайте вращение и как кватернион и как матрицу вращения. Затем преобразуйте кватернион и матрицу вращения к Углам Эйлера в градусах для сравнения.

Задайте силу магнитного поля в микротесла в Бостоне, MA, когда указано истинный север.

m = [19.535 -5.109 47.930];
a = [0 0 9.8];

Определите кватернион и матрицу вращения, которая способна к вращению кадра с магнитного севера на истинный север. Отобразите результаты для сравнения.

q = ecompass(a,m);
quaterionEulerAngles = eulerd(q,'ZYX','frame')

r = ecompass(a,m,'rotmat');
theta = -asin(r(1,3));
psi = atan2(r(2,3)/cos(theta),r(3,3)/cos(theta));
rho = atan2(r(1,2)/cos(theta),r(1,1)/cos(theta));
rotmatEulerAngles = rad2deg([rho,theta,psi])
quaterionEulerAngles =

   14.6563         0         0


rotmatEulerAngles =

   14.6563         0         0

Используйте ecompass, чтобы определить вектор силы тяжести на основе данных из вращения IMU.

Загрузите данные об инерционном модуле измерения (IMU).

load 'rpy_9axis.mat' sensorData Fs

Определите ориентацию корпуса датчика относительно локального кадра NED в зависимости от времени.

orientation = ecompass(sensorData.Acceleration,sensorData.MagneticField);

Чтобы оценить вектор силы тяжести, сначала вращайте показания акселерометра от каркаса кузова датчика до кадра NED с помощью вектора кватерниона orientation.

gravityVectors = rotatepoint(orientation,sensorData.Acceleration);

Определите вектор силы тяжести как в среднем восстановленные векторы силы тяжести в зависимости от времени.

gravityVectorEstimate = mean(gravityVectors,1)
gravityVectorEstimate = 1×3

    0.0000   -0.0000   10.2102

Плавьте смоделированный акселерометр и данные о гироскопе, чтобы отследить вращающуюся платформу с помощью и идеализированных и реалистических данных.

Сгенерируйте траекторию наземной истины

Опишите ориентацию наземной истины платформы в зависимости от времени. Используйте Систему kinematicTrajectory object™, чтобы создать траекторию для платформы, которая не имеет никакого перевода и вращений о его оси z.

duration = 12;
fs = 100;
numSamples = fs*duration;

accelerationBody = zeros(numSamples,3);

angularVelocityBody = zeros(numSamples,3);
zAxisAngularVelocity = [linspace(0,4*pi,4*fs),4*pi*ones(1,4*fs),linspace(4*pi,0,4*fs)]';
angularVelocityBody(:,3) = zAxisAngularVelocity;

trajectory = kinematicTrajectory('SampleRate',fs);

[~,orientationNED,~,accelerationNED,angularVelocityNED] = trajectory(accelerationBody,angularVelocityBody);

Модель, получающая данные IMU

Используйте Системный объект imuSensor, чтобы подражать данным, полученным от IMU, который содержит идеальный магнитометр и идеальный акселерометр.

IMU = imuSensor('accel-mag','SampleRate',fs);
[accelerometerData,magnetometerData] = IMU(accelerationNED, ...
                                           angularVelocityNED, ...
                                           orientationNED);

Плавьте данные IMU, чтобы оценить ориентацию

Передайте данные об акселерометре и данные о магнитометре к функции ecompass, чтобы оценивать ориентацию в зависимости от времени. Преобразуйте ориентацию в Углы Эйлера в градусах и постройте результат.

orientation = ecompass(accelerometerData,magnetometerData);
orientationEuler = eulerd(orientation,'ZYX','frame');

timeVector = (0:numSamples-1).'/fs;

figure(1)
plot(timeVector,orientationEuler)
legend('z-axis','y-axis','x-axis')
xlabel('Time (s)')
ylabel('Rotation (degrees)')
title('Orientation from Ideal IMU')

Повторите эксперимент с реалистической моделью датчика IMU

Измените параметры Системного объекта IMU, чтобы аппроксимировать реалистические данные о датчике IMU. Сбросьте IMU и затем вызовите его с тем же ускорением наземной истины, угловой скоростью и ориентацией. Используйте ecompass, чтобы плавить данные IMU и построить результаты.

IMU.Accelerometer = accelparams( ...
    'MeasurementRange',20, ...
    'Resolution',0.0006, ...
    'ConstantBias',0.5, ...
    'AxesMisalignment',2, ...
    'NoiseDensity',0.004, ...
    'BiasInstability',0.5);
IMU.Magnetometer = magparams( ...
    'MeasurementRange',200, ...
    'Resolution',0.01);
reset(IMU)

[accelerometerData,magnetometerData] = IMU(accelerationNED,angularVelocityNED,orientationNED);

orientation = ecompass(accelerometerData,magnetometerData);
orientationEuler = eulerd(orientation,'ZYX','frame');

figure(2)
plot(timeVector,orientationEuler)
legend('z-axis','y-axis','x-axis')
xlabel('Time (s)')
ylabel('Rotation (degrees)')
title('Orientation from Realistic IMU')

Входные параметры

свернуть все

Показания акселерометра в системе координат корпуса датчика в m/s2, заданном как N-by-3 матрица. Столбцы матрицы соответствуют x - y - и z - оси корпуса датчика. Строки в матрице, N, соответствуют отдельным выборкам. Показания акселерометра нормированы перед использованием в функции.

Типы данных: single | double

Показания магнитометра в системе координат корпуса датчика в µT, заданном как N-by-3 матрица. Столбцы матрицы соответствуют x - y - и z - оси корпуса датчика. Строки в матрице, N, соответствуют отдельным выборкам. Показания магнитометра нормированы перед использованием в функции.

Типы данных: single | double

Формат раньше описывал ориентацию, заданную как 'quaternion' или 'rotmat'.

Типы данных: char | string

Выходные аргументы

свернуть все

Ориентация, которая может вращать количества от глобальной системы координат до системы координат тела, возвратилась как вектор кватернионов или массива. Размер и тип orientation зависят от формата, используемого, чтобы описать ориентацию:

  • 'quaternion'N-by-1 вектор кватернионов с тем же базовым типом данных как вход

  • 'rotmat' – 3 3 N массивом совпадающий тип данных как вход

Типы данных: quaternion | single | double

Алгоритмы

Функция ecompass возвращает кватернион или матрицу вращения, которая может вращать количества от родительского элемента (NED) кадр дочернему элементу (датчик) кадр. Для обоих выходных форматов ориентации оператор вращения определяется путем вычисления матрицы вращения.

Матрица вращения сначала вычисляется с посредником:

R=[(a×m)×aa×ma]

и затем нормированный по столбцам. a и m являются входом accelerometerReading и входом magnetometerReading, соответственно.

Чтобы понять матричное вычисление вращения, рассмотрите произвольную точку на Земле и ее соответствующем локальном кадре NED. Примите каркас кузова датчика, [x, y, z], с тем же источником.

Вспомните, что ориентация корпуса датчика задана как оператор вращения (матрица вращения или кватернион) требуемый вращать количество от родительского элемента (NED) кадр дочернему элементу (корпус датчика) кадр:

[R][pродительский элемент]=[pдочерний элемент]

где

  • R является 3х3 матрицей вращения, которая может быть интерпретирована как ориентация дочернего кадра.

  • Родительский элемент p является вектором 3 на 1 в родительском кадре.

  • Дочерний элемент p является вектором 3 на 1 в дочернем кадре.

Для стабильного корпуса датчика акселерометр возвращает ускорение из-за силы тяжести. Если корпус датчика отлично выравнивается с системой координат NED, все ускорение из-за силы тяжести приезжает z - ось и чтения акселерометра [0 0 1]. Считайте матрицу вращения требуемой вращать количество от системы координат NED до количества, обозначенного акселерометром.

[r11r21r31r12r22r32r13r23r33][001]=[a1a2a3]

Третий столбец матрицы вращения соответствует чтению акселерометра:

[r31r32r33]=[a1a2a3]

Магнитометр, читая точки к магнитному северу и находится в N-D плоскость. Снова, считайте каркас кузова датчика выровненным с системой координат NED.

По определению E - ось перпендикулярна N-D плоскость, поэтому ND = E, в рамках некоторого амплитудного масштабирования. Если каркас кузова датчика выравнивается с NED, и ускоряющий вектор от акселерометра и вектор магнитного поля от магнитометра лежат в N-D плоскость. Поэтому ma = y, снова с некоторым амплитудным масштабированием.

Считайте матрицу вращения требуемой вращать NED к дочернему кадру, [x y z].

[r11r21r31r12r22r32r13r23r33][010]=[a1a2a3]×[m1m2m3]

Второй столбец матрицы вращения соответствует векторному произведению чтения акселерометра и чтения магнитометра:

[r21r22r23]=[a1a2a3]×[m1m2m3]

По определению матрицы вращения, столбец 1 является векторным произведением столбцов 2 и 3:

[r11r12r13]=[r21r22r23]×[r31r32r33]=(a×m)×a

Наконец, матрица вращения нормирована по столбцам:

Rij=Riji=13Rij2,j

Примечание

Алгоритм ecompass использует магнитный север, не истинный север, для системы координат NED.

Ссылки

[1] Fusion датчика с открытым исходным кодом. https://github.com/memsindustrygroup/Open-Source-Sensor-Fusion/tree/master/docs

Расширенные возможности

Генерация кода C/C++
Генерация кода C и C++ с помощью MATLAB® Coder™.

Введенный в R2018b