Обработка информации в системе ориентации на основе МЭМС

УДК 629.7.05

В.Л. Волков, Н.В. Жидкова ОБРАБОТКА ИНФОРМАЦИИ В СИСТЕМЕ ОРИЕНТАЦИИ НА ОСНОВЕ МЭМС

Арзамасский политехнический институт

Предмет, тема, цель работы: Исследование алгоритма бесплатформенной системы ориентации (БСО) с использованием микроэлектромеханических сенсоров (МЭМС).

Методология проведения работы: Составлена математическая модель БСО на основе МЭМС, обеспечивающая вычисление углов ориентации подвижного объекта с использованием кинематических уравнений Пуассона с помощью кватернионов. Использована методика компьютерного эксперимента на основе Matlab для исследования оценок эффективности БСО в условиях дрейфа гироскопов.

Результаты: Исследована работоспособность БСО, использующей МЭМС гироскопы, акселерометры и магнитометры. Приведены результаты по значениям среднего квадрата ошибки (СКО) и предельной ошибки для каждого из углов ориентации при наличии дрейфа гироскопов.

Значения исследования: Применение предложенной методики компенсации дрейфа гироскопов в алгоритме БСО дает приемлемые для практики ошибки.

Новизна/ценность: Применен современный алгоритм БСО, основанный на использовании параметров Родри-га-Гамильтона (кватернионов). Исследование БСО проведено в среде Simulink Matlab, что дает визуальную наглядность схем, возможность их повторения и модификации математических моделей.

Ключевые слова: МЭМС, бесплатформенная система ориентации, кватернионы, средний квадрат ошибки, моделирование, эффективность, Matlab.

Введение

Использование МЭМС в инерциальных измерительных системах сдерживается их относительно невысокими точностными характеристиками [1]. Тем не менее, достижения в МЭМС технологии дают для военной, аэрокосмической техники и техники специального назначения датчики инерции уже с достаточно приемлемыми характеристиками. В последние годы точность инерциальных МЭМС датчиков неуклонно повышается. Нестабильность смещения гироскопов компании Silicon Sensing теперь составляет не более 0,3 °/ч, а акселерометра - 50 м^ [2]. Более значительный эффект для инерциальных измерений дает применение инерциальных блоков, т.е. совмещение в одном корпусе 3-осевых датчиков угловых скоростей (ДУС) с акселерометрами, магнитометрами и датчиками температуры. В этом случае путем совместной калибровки их характеристик удается достичь высокой точности измерения инерциальных параметров. Однако следует отметить недоступность зарубежных электронных компонентов, к которым относятся и МЭМС, для применения в военной технике и технике специального назначения, что заставляет искать альтернативные - алгоритмические пути решения проблем точности инерциальных навигационных систем.

Особенности построения бесплатформенных инерциальных навигационных систем (БИНС) и их составных частей - бесплатформенных систем ориентации (БСО) стали особенно актуальны в условиях интенсивно развивающейся техники мобильных подвижных устройств. Системы ориентации стали необходимой составной частью беспилотных летательных аппаратов, роботов, мобильных навигационных устройств, подвижных объектов специального назначения. Бесплатформенный принцип исполнения БСО стал практически основным при использовании МЭМС и микро ЭВМ в мобильных подвижных объектах.

В последние годы появилось множество работ в области бесплатформенных инерци-альных систем, это связано, в первую очередь с интенсивным развитием летательных аппаратов 4-го и 5-го поколений, беспилотных летательных аппаратов и робототехники. Здесь следует отметить фундаментальную работу В.Н. Бранца и И.П. Шмыглевского [3]. Осново-

© Волков В.Л., Жидкова Н.В., 2015.

полагающими также являются работы [4-6]. Особенно следует отметить современное учебное пособие В.Я. Распопова и В.В. Матвеева [4], эффективно используемое для изучения основ БИНС.

Математическое представление БСО, как правило, реализуется в компьютерные модели с целью проверки работоспособности алгоритмов БСО и с целью исследования их эффективности еще на этапе проектирования. Компьютерное исследование алгоритмов БСО наиболее удачно реализуется на основе инженерного программного обеспечения Ыа^аЪ БтиМпк [4, 7].

Основная задача данной статьи состоит в исследовании БСО, реализуемой на основе параметров Родрига-Гамильтона (кватернионов) путем компьютерного моделирования в Ыа^аЪ 81шиИпк. Следует отметить, что большинство авторов по БИНС и БСО в настоящее время отдают предпочтение именно алгоритмам на базе математических моделей с параметрами Родрига-Гамильтона [3 - 6]. Хорошей особенностью БСО на основе кватернионов является линейность и интегрируемость кинематических уравнений при любых углах Эйлера-Крылова.

Базовой основой для данной статьи явились исследовательские материалы, представленные в [6]. В данной статье рассматривается составление схем БСО в инструментальной среде 8тиМпк с целью проверки математических моделей и получения результатов исследования эффективности БСО при их работе в различных динамических режимах. Исполнение математических моделей БСО в 8тиМпк позволяет получить более понятную визуальную модель системы и реализовать многочисленные операции дифференцирования и интегрирования более гибко на стандартных элементах. Важное место в статье уделено методике и эффективности компенсации основной погрешности БСО, возникающей вследствие дрейфа гироскопов.

Статья имеет, прежде всего, практическую направленность, позволяющую повторить предлагаемые модели и схемы, получить результаты симуляции и оценить эффективность алгоритмов БСО. Работу над статьей авторы проводили на основе многочисленных источников информации, отраженных в списке литературы, и собственных разработок [6].

Математическая модель БСО

Источником информации алгоритма ориентации с кватернионами являются проекции угловой скорости, получаемые на основе сигналов трех датчиков угловых скоростей (ДУС) и вычисляемых проекций абсолютной угловой скорости географического трехгранника. Выходом алгоритма ориентации являются оценки углов тангажа, крена и рыскания: и, у, у.

Алгоритм определения ориентации по данным об угловой скорости объекта основан на

решении двойного кинематического уравнения Пуассона [4]:

•

2 Л = Л о о Л + Л

где Л - кватернион, определяющий переход из географической системы координат OXgYgZg к связанной системе ОХУ^, й^, - гиперкомплексные отображения векторов абсолютной угловой скорости в связанной и горизонтальной системах координат соответственно;

- норма кватерниона; Х0, Х1, Х2, Х3 - параметры Родрига-Гамильтона.

Абсолютная угловая скорость ю£=[юхз,юге,ю2$], измеряемая ДУСами в связанной системе координат, складывается из вектора переносной угловой скорости земной системы координат, вызванной вращением Земли, из вектора переносной угловой скорости, связанной с перемещением объекта относительно Земли с линейной скоростью и из вектора угловой скорости связанной системы координат относительно базовой [4]. В векторно-матричном виде это соответствует соотношению:

ю^ = Люа + ю0,

где юа=[юх^ ю^, ю^] - проекции переносной угловой скорости в земной системе координат;

ю0=[юх, юг, ю^ - проекции относительной угловой скорости в связанной системе координат; A - матрица направляющих косинусов, имеющая вид:

A =

— i

cos \ cos и sin и - sin \ cos и

cos \ sin и cos Y + sin \ sin у cos и cos у sin \ sin и cos у + cos \ sin у .

cos \ sin и sin у + sin \ cos у - cos и sin у - sin \ sin и sin у + cos \ cos у

Относительная угловая скорость ю0 определяется следующим соотношением [4]:

= \\ sin и + у;

ш7 = \\ cos и cos у + и sin у;

ш2 = cos и sin у + и cos у; где и, у, у - углы Эйлера-Крылова: тангаж, крен и рыскание; у\ ^ - угловые скорости соответственно по тангажу, крену и курсу.

Проекции Юха, yg, zg, в свою очередь, вычисляются по следующим известным соотношениям [4]:

юж = U cos ф + Vza /R;

7G = Usin ф + tanф-^j /R;

zg =-VXG / R-

где Vzg, Vxg - составляющие линейной скорости объекта соответственно вдоль оси Zg и Xg географической системы координат; U - угловая скорость вращения Земли; ф - широта местоположения объекта; R - радиус Земли.

Решение уравнения (1) позволяет определить собственный кватернион ориентации Л между горизонтальной и связанной системами координат при задании начальных условий Л(0) [4]:

(0) = cos(\ /2)cos(u0 /2)cos(^ /2) - sin(\0 /2)sin(u0 /2)sin(^ /2);

(0) = cos(\0 / 2) cos(u0 / 2) sin(у0 /2) + sin(\0 / 2) sin(u0 / 2) cos^0 / 2); ^

Я2(0) = sin(\0 /2)cos(u0 /2) cos^0 /2) + cos(\0 /2)sin(u0 /2)sin^0 /2);

(0) = cos(\ / 2) sin(u0 / 2) cos(^ /2) - sin(\0 / 2) cos(u0 / 2) sin(у0 / 2), где у0, u0, y0 - начальные значения углов рыскания, тангажа и крена.

Начальные значения углов и0 и у0 определяются на основе проекций ускорения свободного падения на оси связанной системы координат (требуется информация с акселерометров):

и0 = arcsin(-gx / g); у0 = arctan(gz / g7 ),

где gX, gY, gz - проекции ускорения свободного падения на оси связанной системы координат, измеряемые акселерометрами при равномерном движении объекта.

Для вычисления угла требуется информация с трехосного магнитометра. Магнитометр измеряет составляющие напряженности магнитного поля Земли по осям связанной системы координат: HX, HY, HZ. Значение H для магнитного поля Земли составляет: H=0,25 - 0,65 Oe (эрстед).

В качестве трехосевых магнитометров часто применяют магниторезистивные датчики фирмы Honeywell [8]. Для построения электронного компаса Honeywell выпускает набор HMC1055, включающий двух- и одноосевой датчики магнитного поля, а также М-датчик ускорения для определения крена объекта. HMC1055 является фактически трехосевым с диапазоном ±6 Гаусс и разрешающей способностью 85 мкГаусс. Вычисление угла курса в диапазоне почти ±90 град. на основе показаний магнитометра осуществляется на основе соотношения:

= arcsin{

где HY, HZ - составляющие напряженности магнитного поля Земли по осям связанной системы координат; D - магнитное склонение.

Вычисление начальных условий для кинематического уравнения (1) необходимо производить периодически для компенсации довольно значительного дрейфа показаний ДУСов (0,02 - 0,03 град/с). Однако режим компенсации можно применять только в моменты прямолинейного равномерного движения объекта. Такой режим можно определить по показаниям

акселерометров, когда ^a2х + a2y + a2 « g, где ax2, ay2, a 2 - составляющие кажущегося ускорения на оси связанной системы координат.

По элементам найденного в результате интегрирования кватерниона Л (см. уравнение (1)) на конечном этапе определяются углы ориентации [4]:

и = arcsin(2À1À 2 + 2Я0Я3);

у = arctan

^ = arctan

Компенсация дрейфа показаний ДУСов осуществляется на основе показаний акселерометров, измеряющих кажущееся ускорение в виде проекций на оси связанной системы координат и показаний магнитометров. Известные алгоритмы выделения ускорения свободного падения из кажущегося ускорения и информация о скорости движения объекта относительно земной системы координат позволяют дополнительно определить два угла поворота связанной системы координат относительно земной (см. соотношения (7)) и еще один угол (курс) определяется по показаниям трехосного магнитометра (8). Далее применяется компенсационная схема на основе ПИ-регулятора [9], которая позволяет отфильтровывать имеющиеся измерительные шумы акселерометров и выделить и компенсировать ошибки дрейфа гироскопов.

Схема для симуляции БСО составлена в инструментальной среде Matlab Simulink [10] и представлена на рис. 1. Она имеет блоки: Objects Model, Wg, AnglO, Poisson, Quaternion to Angle, Compensator, Result и ряд дополнительных.

Рис. 1. Схема симуляции БСО на основе кватернионов

При моделировании алгоритма БСО в Simulink необходима имитация угловых движений подвижного объекта и показаний сенсоров, измеряющих это движение. Для имитации истинных значений углов и, у, у, угловых скоростей roX, roY, roZ и вектора проекций ускорения свободного падения gX, gY, gZ создан блок Objects Model. Этот блок вырабатывает векторы Angl=[u,y,y], W=[roX,My,Mz], G=[gx, gY, gZ] и вектор H=[Hx,Hy,Hz].

Схема блока Objects Model представлена на рис. 2. Истинные угловые движения объекта имитируются источниками с названиями psi, tet, gam. В данном примере моделирования приняты законы изменения углов рыскания, тангажа и крена синусоидальными со следующими параметрами: рыскание - начальное значение yo =arctg(Vzg/Vxg), амплитуда Av=1°, ча-

tet Angle Velosty

gam h

tet gam W -psi

t— psi —

Acselerometer

стота wy=0,1 Гц; тангаж - начальное значение uo=n/24, амплитуда Аи=1,5о, частота w^=0,15 Гц; крен - начальное значение yo=rc/24, амплитуда Ау=0,3о, частота wY=0,2 Гц.

Выработка вектора проекций ускорения свободного падения G=[gx,gY,gz] в схеме рис. 2 производится субсистемой - Accelerometer на основе формул в соответствии с работой [5]: gx=gsin(v), gY=gcos(y) cos(u), gz= gsin(y)cos(u). Для имитации идеальных проекций относительной угловой скорости на оси связанной системы координат: roX, roY, roz на основе истинных и, у, у создан блок Angle Velosity, работающий на основе уравнений (4). Блок Mag осуществляет имитацию показаний магнитометра на основе формул из [5]:

Нх = L cos у cos и + B sin и;

HY = L sin у sin у + B(cos и - L cos у sin и) cos у;

Н = L sin у cos у - B(cos и - L cos у sin и) sin у.

W

OD СЮ

G

сю

P

Рис. 2. Блок Objects Model

где ,8=Hsin(I); L=Hcos(I) - вертикальная и горизонтальная составляющие напряженности маг-

2 2 2 1/2

нитного поля Земли; #=(#х +#y +#z ) ' - полная напряженность магнитного поля Земли.

Входными данными для блока Mag являются истинные значения углов Эйлера-Крылова Angl. Блок Mag выдает показания магнитометров в виде вектора H. Далее в блоке Psi эти показания используются для вычисления угла курса в соответствии с формулой (8). Угол курса у0 является выходом блока Psi. Выработка угла магнитного наклонения I осуществляется на основе информации о широте и долготе местности по карте магнитных наклонений. В данном случае имитатора показаний магнитометра, используя истинные значения углов и, у, у, при вычислении угла курса в соответствии с формулой (8) получаем угол по отношению к географическому северу, т.е. то, что используется в навигации и ориентации объекта в земной системе координат. Для определенности при имитации H выбрано конкретное значение магнитной напряженности Земли H=0,0314 [А/м] (в системе СИ -в амперах на метр), соответствующее местоположению объекта.

Блок имитации показаний ДУСов Measure содержит схему вычисления абсолютной угловой скорости, измеряемой ДУСами в соответствии с матричным соотношением (2). Его входами являются векторы ro0=[roX,roY,roZ], Ang=[u,y,y], <c=[<xg,<Yg,<Zg] и вектор ошибок dW=[droX,droY,droZ]. Этот блок вырабатывает вектор измерений roS=[roXS,roYS,roZS]. Для каждого ДУС добавлен соответствующий дрейф: dwx, dwy, dwz: WS=<S+dW. Блок Wg формирует вектор переносной угловой скорости в земной системе координат <G=[<xg, <Yg, <Zg].

Для реализации рассматриваемого алгоритма БСО на кватернионах в Simulink требуется ряд матричных блоков - один из них производит умножение кватернионов Л и Q в уравнении Пуассона. Для формирования матрицы M(X) используется стандартный блок Matrix Concatenation, объединяющий столбцы в матрицу. Входным параметром блока является кватернион L, а выходным квадратная гиперкватернионная матрица M(X). Блок формирования матрицы гиперкватерниона типа M(X) назван Matrix Quaternion.

Умножение кватернионов Л и Q в уравнении Пуассона соответствует умножению следующих матриц [4]

X 0 -Xi -X 2 -X3" " 0 "

M (X) = Xl X 2 X 0 X 3 -X3 X0 X 2 -Xi ; ra = ю7

_Xs -X 2 Xl X 0 _ _

Для формирования вектора кватерниона угловой скорости на основе второй формулы (11) используется простая схема, добавляющая ноль в верней части трехмерного вектора угловой скорости. Перемножение матриц Ы(Х) и Q выполняется с помощью стандартного блока Matrix Multiply.

Для составления схемы решения уравнения Пуассона (1) требуется еще интегратор, несколько масштабирующих элементов и несколько сумматоров. Для задания начальных условий интегрирования потребуется еще блок QuatO, вычисляющий кватернион ¿(0). Блок решения уравнения Пуассона (1) на основе описанных элементов представлен на рис. 3.

Рис. 3. Блок решения уравнения Пуассона (1) - Poisson

Этот блок назван Poisson. Входами этого блока являются трехмерные векторы угловых скоростей W и Wg, соответственно связанной и географической систем координат. Выходом блока является кватернион, являющийся решением уравнения Пуассона (1).

В верхней части схемы рис. 3, справа представлена схема формирования элемента коррекции нормы кватерниона Л(1-Л ). В итоге, при работе схемы рис. 3 на ее выходе получается кватернион оценки вектора углов Эйлера-Крылова. Его уже можно преобразовать в углы по формулам (9). Блок QuatO, входящий в схему рис. 3 и вычисляющий L(0) - элементы кватерниона в начальный момент времени, работает на основе формул (6). Входами этого блока являются начальные значения углов тангажа, крена и рыскания - uo, yo, yo, а выходом кватернион начальных условий, элементы которого вычисляются на основе формул (6).

Блок, преобразующий кватернион L в углы по формулам (9), на схеме рис. 1 назван Quaternion to Angle. Блок Compensator предназначен для компенсации ошибок дрейфа гироскопов. В основе блока Compensator лежит компенсационная схема на основе ПИ-регулятора [9]. Ошибки оценивания углов Эйлера относительно дополнительных углов Ao, вычисленных на основе данных акселерометров и магнитометров, подаются на ПИ-регуляторы. При этом выделяются низкочастотные ошибки дрейфа гироскопов. Выделенные ошибки вычитаются из неточных оценок, что дает существенно более точные результаты оценивания углов Эйлера.

Блок Result схемы рис. 1 предназначен для регистрации результатов симуляции БСО.

Результаты симуляции БСО с кватернионами

Выполняя симуляцию без подключения блока компенсации дрейфа гироскопов, построим также графики оценивания углов Эйлера на фоне реальных углов в случае с дрейфом гироскопов. Для примера дрейф гироскопов задан следующим: dwx=0,01 град/с; dwy=0,01 град/с; dwz=0,01 град/с. Графики оценок углов Эйлера на фоне их истинных значений представлены на рис. 4.

Средние квадраты ошибок, рассчитанные по результатам эксперимента, получены следующие: СКОт=0,23 - для угла тангажа, СКОР=0,1 - для угла рыскания, СКОа=0,17 - для угла крена. Пересчет в предельные погрешности (при Av=1°, Au=1,5°, AY=0,3°) дает следующие результаты в процентах: dT%=1834%, dP%=1819%, dG%=7934%.

Рис. 4. Графики оценок углов Эйлера на фоне их истинных значений

При включенном блоке компенсации дрейфа гироскопов и достаточно большом уровне дрейфа (^х=0,01 град/с; ^ну=0,01 град/с; ^¿=0,01 град/с) получены графики оценок углов Эйлера на фоне их истинных значений (см. рис. 5) и численные результаты погрешностей.

0.03 0.02 0.01 0

-0.01 -0.02 -0.03

0

50

100 40

50

100

0.57 0.56 0.55 0.54 0.53 0.52

0

50

100

Рис. 5. Графики оценок углов Эйлера на фоне их истинных значений с учетом компенсации дрейфа гироскопов

Средние квадраты ошибок, рассчитанные по результатам эксперимента, в этом случае

12 13

получены следующие: СК0т=1,710- - для угла тангажа, СК0Р=6,5 10- - для угла рыскания, СКОа=7,2 10"13 - для угла крена. Пересчет в предельные погрешности (при Ау=1о, Аи=1,5о, Ау=0,3°) дает следующие результаты в процентах: ^7%=0,005%, оР%=0,0046%, ¿а%=0,015%.

Выводы

Математическая модель бесплатформенной системы ориентации построена на основе традиционных кинематических уравнений Пуассона с использованием параметров Родрига-Гамильтона. Применение алгоритма БСО на кватернионах имеет существенные преимущества по сравнению с другими алгоритмами (на углах Эйлера-Крылова, на направляющих косинусах), так как имеет существенно меньшие погрешности расчетов при прочих равных условиях.

Моделирование алгоритма БСО проведено с помощью программной среды МайаЪ иНпк. В качестве результатов моделирования приведены графики основных процессов и численные результаты по ошибкам оценивания углов ориентации. Результаты компьютерного эксперимента подтверждают эффективность алгоритма БСО в системах определения пространственного положения маневренных подвижных объектов в широком диапазоне сигналов и шумов. Ввиду незначительного расхождения результаты в виде графиков оценок практически совпадают с реальными углами, что доказывает эффективность компенсации ошибок дрейфа показаний гироскопов. В итоге, компенсация дрейфа гироскопов, применяемая в алгоритме оценивания на основе кватернионов дает высокую точность оценок (не хуже 0,02%).

Новым результатом, полученным в статье, является алгоритм компенсации дрейфа в измерительных данных гироскопов. Эффект компенсации дрейфа достигнут с помощью информация с микроакселерометров и магнитометров при использовании схемы с отрицательной обратной связью на основе ПИ-регуляторов. Таким образом, решена задача компенсации дрейфа показаний микрогироскопов, как основной ошибки БСО.

Библиографический список

1. Сысоева, С.С. Тенденции рынка Higt-end МЭМС датчиков инерции. Новые уровни характеристик и исполнения // Компоненты и технологии. 2014. № 4. С. 41-46.

2. Бакмачев, А. МЭМС-гироскопы и акселерометры Silicon Sensing: английские традиции, японские технологии // Компоненты и технологии. 2014. № 4. С. 18-26.

3. Бранец, В. Н. Применение кватернионов в задачах ориентации твердого тела / В. Н. Бранец, И. П. Шмыглевский. - М.: Наука, 1973. - 320 с.

4. Матвеев, В. В. Основы построения бесплатформенных инерциальных навигационных систем. / В. В. Матвеев, В. Я. Распопов; под ред. В. Я. Распопова. - СПб.: ГНЦ РФ ОАО «Концерн «ЦНИИ «Электроприбор», 2009. - 280 с.

5. Распопов В.Я. Микросистемная авионика: учеб. пособие / В.Я. Распопов. - Тула: «Гриф и К», 2010. - 248 с.

6. Жидкова, Н.В. Моделирование бесплатформенной системы ориентации / Н.В. Жидкова, В.Л. Волков // Современные проблемы науки и образования. 2015. № 1; URL http: //www .science -education.ru/121-17099.

7. Белочкин, П.Е. Особенности моделирования бесплатформенной системы ориентации по уравнениям Эйлера в среде Matlab и Mathcad / П. Е. Белочкин , Д. А. Кацай // Современные наукоемкие технологии. 2014. № 5. С. 18-20;

8. Маргелов, А. Магниторезистивные датчики положения компании Honeywell. Chip News. URL: www.chip-news.ru. 2005. № 95. С. 52-56.

9. Premerlani, W. Direction Cosine Matrix IMU: Theory / W. Premerlani, P. Bizard. URL: http://gentlenav.googlecode.com/files/DCMDraft2.pdf. 2009. 30 p.

10.Дьяконов, В.П.Matlab. Полный самоучитель / В.П. Дьяконов. - М.: ДМК Пресс, 2011. -768 с.

Дата поступления в редакцию 13.06.2015

V.L. Volkov, N.V. Zhidkova

INFORMATION PROCESSING IN SYSTEM OF ORIENTATION,

BASED ON MEMS

Arzamasskiy pollytechnic institute (branch), Nizhny Novgorod state technical university n.a. R.E. Alexeev

Purpose: Study algorithm strapdown orientation system (SSR) using micro electro-mechanical sensors (MEMS). Design/methodology/approach: The mathematical model SSR MEMS calculates orientation angles of moving object using Poisson kinematic equations using quaternions. Used computer simulation method based on Matlab for studying effectiveness evaluations SSR in a gyro drift.

Results: A study using SSR performance MEMS gyroscopes, accelerometers and magnetometers. The results of values of mean square error (RMSE) and limit of error for each of angles of orientation in presence of drift gyros. Research limitations/implications: Application of proposed method gyroscopes drift compensation algorithm gives BSO acceptable for practice of error.

Originaly/value: Using modern algorithm SSR, based on of parameters Rodrigues-Hamilton (quaternions). A study conducted among SSR Simulink Matlab, giving visual clarity schemes, possibility of repetition and modification of mathematical models.

Key words: MEMS strapdown orientation system, quaternions, the mean square error, modeling, efficiency,

Matlab.