Использование одноосного гироскопа для определения ориентации макета в лабораторных условиях

( Determination of the Orientation of a Suspended Mock-up Using One-Axis Gyroscope,
Preprint, Inst. Appl. Math., the Russian Academy of Science)

Иванов Д.С., Овчинников М.Ю.
(D.S.Ivanov, M.Yu.Ovchinnikov)

ИПМ им. М.В.Келдыша РАН

Москва, 2008
Работа выполнена при частичной поддержке Российского фонда фундаментальных исследований (проект № 06-01-00389) и DAAD (программа Leonhard Euler, реферат 325)

Аннотация

Рассматривается задача определения ориентации твердого тела, подвешенного на струне, с использованием одноосного датчика угловой скорости (гироскопа). Изложены элементы теории расширенного адаптивного фильтра Калмана. Дано описание фильтров Калмана, основанных на различных моделях движения тела. Проводится компьютерное моделирование процесса определения ориентации тела с помощью фильтров Калмана. Дано описание макета, созданного в ИПМ им. М.В.Келдыша РАН, на котором реализована работа фильтров Калмана, использующих в качестве входных данных показания гироскопа.

Abstract

A problem of attitude determination of a string-suspended muck-up equipped with one-axis angular velocity sensor (gyroscope) is considered. Particulars of the extended adaptive Kalman filter theory are stated. Kalman filters based on various models of the mock-up are presented. Computer simulation of the attitude determination of the mock-up is carried out using the filters. Laboratory facility for laboratory simulation is developed at the KIAM RAS which is used for the filter verification.

Содержание

1. Введение                                                                                          3

1.1 Способы определения ориентации тела                                    3

1.1.1 Позиционные датчики                                                      3

1.1.2 Датчики угловых скоростей                                            5

1.2 Методы устранения шума датчика                                             6

1.2.1 Статистические методы                                                   6

1.2.2 Фильтрация                                                                      7

2. Элементы теории фильтра Калмана                                          7

2.1 Линейная задача. Основные понятия                                        7

2.2 Коэффициент обратной связи                                                    9

2.3 Коррекция значения ковариационной матрицы ошибки        10

2.4 Фильтр Калмана для нелинейных систем                                 10

2.5 Использование фильтра Калмана                                              11

2.6 Адаптивный фильтр. Метод “Уточнение ковариации”            12

3. Уравнения движения макета                                                      13

3.1 Анализ сил и моментов, действующих на макет                       13

3.2 Вывод уравнений движения тела на струне                              14

4. Применение фильтра Калмана                                                    16

4.1 Перманентное вращение                                                            16

4.2 Крутильные колебания                                                              17

5. Математическое моделирование определения ориентации   18

5.1 Структура программы                                                              18

5.2 Результаты моделирования                                                       20

6. Проведение измерений с использованием макета                    23

6.1 Схема макета                                                                              24

6.2 Характеристика и устройство датчиков угловой скорости     25

6.3 Программное обеспечение                                                        26

6.3.1 Работа программы                                                          26

6.3.2 Интерфейс                                                                        27

6.4 Результаты экспериментов                                                        28

7. Заключение                                                                                    30

8. Благодарности                                                                              30

9. Список литературы                                                                      31



1. Введение

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

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

1.1Способы определения ориентации тела

Измерительные средства, с помощью которых возможно определить ориентацию тела, можно разделить на две группы датчиков, которые измеряют непосредственно или косвенно (через пересчет) одну из компонент вектора фазового состояния. Первая группа – это датчики направления, определяющие в связанной системе некоторое известное в неподвижной системе координат направление. Вторая группа – датчики угловой скорости, измеряющие непосредственно вектор угловой скорости тела. Рассмотрим подробнее, как с помощью конкретных измерительных средств определяют ориентацию тела.

1.1.1 Позиционные датчики

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

·                 Звёздная камера

Подобные устройства, используя привязку наблюдаемых звёзд к известному звёздному каталогу, обеспечивают определение трёхосной ориентации аппарата относительно инерциального пространства. Этот датчик представляет собой видеокамеру, снимающую звёздное небо. В рассматриваемой в [1] макетной реализации камеры результатом съёмки является видеофайл. Видеофайл далее преобразуется в последовательность графических файлов. Графические файлы являются входными для алгоритма определения углового положения спутника, версия которого описана в [1]. Основные этапы работы этого алгоритма можно описать следующим образом.

– Определяется положение центров масс всех наблюдаемых звёзд относительно системы координат, связанной с сенсором камеры.

– Определяется угловое положение распознанных звёзд относительно системы координат, связанной с камерой.

– Вычисляются угловые расстояния между наблюдаемыми звёздами.

– В соответствии со звёздным каталогом проводится идентификация наблюдаемых звёзд.

– Определяется направление оптической оси камеры относительно инерциальной и орбитальной системами координат.

– Определяется ориентация спутника относительно инерциальной и орбитальной системами координат.

·                 Солнечный датчик

Основная задача солнечного датчика – это определение направления на Солнце (азимута и угла места по отношению к спутнику). Рассмотрим его устройство на примере солнечного датчика, базирующегося на радиационно-устойчивом APS детекторе [2], закрытом непрозрачной маской с одним (или несколькими равноудаленными) отверстиями диаметром 0.2 мм, которая существенно ограничивает поступающий на детектор поток солнечного излучения. Схема работы датчика представлена на рис.1.1.

Солнечный свет, проходя через отверстие в маске, формирует на плоскости фотодетектора некоторое изображение источника. Положение центра масс изображения дает возможность определить ориентацию прибора относительно направления на источник. Углы ориентации сенсора a и b могут быть определены из простейших геометрических соображений, согласно формулам

Здесь  - положение центра масс изображения в системе координат OXY, связанной с плоскостью фотодетектора, F – расстояние от маски до плоскости фотодетектора,  – координаты центра масс изображения в случае, когда источник находится в зените (фактически – это расчетные значения проекции центра отверстия маски на плоскость фотодетектора).


·                 Магнитометр

Как правило, подобные устройства состоят из трёх индукционных сенсоров, измеряющих окружающее магнитное поле в трех взаимно перпендикулярных направлениях и, тем самым, определяют ориентацию и величину вектора напряженности поля. На околоземной орбите направление магнитного поля Земли известно в каждой точке пространства. В первом приближении геомагнитное поле подобно полю диполя (рис.1.2). Таким образом, зная ориентацию вектора напряженности поля в связанной системе координат, можно определить ориентацию спутника в орбитальной системе координат [3].

С помощью позиционных датчиков, определяющих непосредственно ориентацию тела, можно также оценить и скорость вращения тела, используя различные алгоритмы вычисления скорости изменения определяемого направления. Однако во многих случаях эти методы не дают удовлетворительной точности. Следует также отметить в числе недостатков звёздной камеры и солнечного датчика, что их можно использовать только при определённых условиях. Так, солнечный датчик не может работать, когда он ориентирован против Солнца или находится в тени Земли, а звёздный датчик не может определять ориентацию, когда он направлен на Землю (не видит звёздного неба) и на Солнце (происходит засветка звёзд). В этих случаях невозможно определить ориентацию спутника.

1.1.2 Датчики угловой скорости

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

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

– Задаётся положение тела относительно некоторой неподвижной системы координат в начальный момент времени.

– Интегрируется кинематическое уравнение

,                                                                        (1.1)

где  - вектор угловой скорости вращения тела, определяемый с помощью датчиков, a  - кватернион перехода из неподвижной системы координат в систему координат, жестко связанную с телом [4]. Таким образом, решается задача определения ориентации.

Точность решения уравнения (1.1) в общем случае зависит от метода реализации этого решения и от ошибок первичной информации, то есть ошибки в определении .

Рассмотрим далее такую схему реализации решения кинематического уравнения, в которой аналоговый сигнал, пропорциональный измеряемой угловой скорости, преобразуется в цифровой и далее интегрирование уравнений выполняется на ЦВМ. Одной из составляющих ошибки будет ошибка из-за дискретности измерений при преобразовании первичной информации. Если квант преобразования равен ε (ε=/ , где N-разрядность преобразователя), то оценка ошибки  дает представление о точности реализации данного метода [5]. Однако, как показывает практика, если разрядность преобразователя достаточно велика (к примеру, для 2-х байтового значения угловой скорости N=16), ошибка реализации метода решения (1.1) очень мала по сравнению с шумом датчика. А неточность в значении угловой скорости после интегрирования уравнения (1.1) даёт накапливающуюся ошибку в определении ориентации.

Следует отметить, что с помощью измерений датчиков угловой скорости в каждый момент времени можно определить полный вектор фазового состояния. При этом датчики угловой скорости не имеют ограничений относительно взаимного положения с Солнцем или Землёй, как солнечный или звёздный датчики. Но у датчиков угловой скорости существует ограничение относительно пределов измерения: датчик не может измерить скорость выше некоторой заявленной, и измерить точнее, чем это позволяет уровень шума. В числе недостатков нужно назвать то, что, так как неизбежна некоторая ошибка в задании положения в начальный момент времени, то дальнейшее определение ориентации происходит с ошибкой.

1.2 Методы устранения шума

Все реальные датчики положения и угловой скорости, выдают не точное значение положения или скорости соответственно, а с некоторой ошибкой. Возникает задача устранения шума датчика – по зашумлённым данным оценить вектор состояния тела, то есть восстановить его реальное движение. Рассмотрим методы устранения шума.

1.2.1 Статистические методы

Статистические методы основаны на минимизации функционала – суммы квадратов разности между вычисленными на основе математической модели движения и измеренными датчиками параметрами. К примеру, на спутнике имеется магнитометр и солнечный датчик. Пусть для известной орбиты и для тех моментов времени, для которых имеются измерения, рассчитывается компоненты векторов направления на Солнце S и направление вектора магнитного поля Н в инерциальной системе отсчёта. В результате получаем набор

            (n=1,…,N; i=1,2,3)

относящихся к моменту времени  результатов измерений величин , в связанной системе координат и расчетные значения , в инерциальной [6]. Введём вектор параметров α, которые необходимо определить, состоящий из начальных значений трех углов Эйлера и трех компонент угловой скорости, а z(t, z0(t0)) – решение уравнений движения спутника относительно центра масс. Это решение получается их численным интегрированием и для него рассчитываются S(z) и Н(z) и матрица  перехода от инерциальной к связанной системе координат. В свою очередь решение уравнений движения z зависит от вектора состояния α, таким образом, z=z(α). Построим функционал

Ф=ФSH,

где

, ,

N – число измерений. Таким образом, оценка вектора α сводится к минимизации функционала Ф(α) по элементам α.

Описанный метод называется методом наименьших квадратов (МНК). МНК требует довольно большого объёма вычислений, и, как следствие, большого времени для расчетов, и поэтому редко применяется на бортовом компьютере спутников для управления движением, а используется для наземной постэкспериментальной обработки информации. Подчеркнём, что МНК плохо работает в случае изменения движения тела под действием каких-либо факторов (к примеру, уменьшилась высота орбиты спутника и на него стала действовать большая аэродинамическая тормозящая сила). В этом случае математическая модель, используемая в МНК, неудовлетворительно описывает движение и в результате оценка вектора α становится неудовлетворительной и нужно менять математическую модель или модифицировать ее параметры.

1.2.2 Фильтрация

Зачастую возникают задачи определения динамической оценки вектора состояния α в условиях шума измерений и шума системы (воздействие на систему различных возмущающих факторов). Такие задачи называются задачами фильтрации.

В настоящей работе рассматривается фильтр Калмана – последовательный рекурсивный алгоритм оценки, использующий принятую модель динамической системы для получения оценки. Преимущество этого алгоритма состоит в том, что он не требует больших вычислительных мощностей, и поэтому может быть реализован на бортовом компьютере спутника. Причем с помощью него можно получать оценку вектора состояния в режиме реального времени, что необходимо, к примеру, для реализации активного управления. Кроме того, фильтр способен подстраиваться в случае неточно выбранных параметров математической модели, на которой он основывается. Этот механизм подробно разобран ниже в разделе 2.6 настоящей работы.

2. Элементы теории фильтра Калмана

Фильтр Калмана – последовательный рекурсивный алгоритм, использующий принятую модель динамической системы для получения оценки, которая может быть существенно скорректирована в результате анализа каждой новой выборки измерений во временной последовательности. Этот алгоритм находит применение в процессе управления многими сложными динамическими системами, например, непрерывными производственными процессами, самолетами, кораблями и космическими аппаратами. При управлении динамической системой, прежде всего, необходимо полностью знать ее фазовое состояние в каждый момент времени. Но измерение всех переменных, которыми необходимо управлять, не всегда возможно и в этих случаях фильтр Калмана является тем средством, которое позволяет восстановить недостающую информацию посредством имеющихся неточных (зашумленных) измерений. При изложении элементов теории фильтра Калмана будем следовать обзорной работе [7].

2.1 Линейная задача. Основные понятия

Предполагаем, что стохастическая система может быть описана моделями динамики и измерений (для непрерывного и дискретного случаев), представленными в табл.2.1

Табл.2.1

Модель

Непрерывное время

Дискретное время

Система

Измерения

Шум

системы

Шум

измерений

Здесь  вектор состояния динамической системы, который является случайным Гауссовским процессом,  - измерения, полученные в момент времени . Шум системы и шум измерений  и  также являются Гауссовскими случайными процессами с нулевым математическим ожиданием. Другими словами мы предполагаем наличие в системе и измерениях белого шума.

Задача фильтрации состоит в том, чтобы найти оценку вектора состояния системы , которую мы будем обозначать , являющуюся функцией измерений и которая минимизирует среднеквадратичную ошибку

,

где M симметричная положительно-определенная матрица.

Подпись: КоррекцияФильтр Калмана работает по системе прогноз-коррекция. Допустим, что в момент времени  получена оценка вектора состояния системы  и теперь хотим получить оценку в момент . Для этого строим прогноз оценки , базируясь на , получаем измерения  и далее корректируем оценку в момент , базируясь на прогнозе и измерениях, и получаем окончательную оценку вектора состояния  (рис. 2.1).  называется априори оценка, называется апостериори оценка.



Рис.2.1. Принцип работы фильтра Калмана

Второй момент случайного процесса может быть описан в терминах ковариационной матрицы

,                                              (2.1)

Матрица P(t) называется ковариационной матрицей ошибки оценки вектора состояния (далее – ковариационная матрица ошибки).

Очевидно, что для получения прогнозируемой оценки  необходимо интегрировать модельное динамическое уравнение  с начальным условием . А каким образом мы можем получить прогнозируемую оценку ковариационной матрицы ошибки? Для того чтобы ответить на этот вопрос рассмотрим непрерывную модель системы

,

приведенную в табл.2.1. Решение этого уравнения с начальным условием x(t0) и переходной матрицей состояния Ф(t, t0) записывается в виде

.

Математическое ожидание этой величины –

.

Тогда

.     (2.2)

Подставляя (2.2) в (2.1), производя некоторые преобразования и вычисляя первую производную функции P(t), в результате получим

.                                                            (2.3)

Это матричное дифференциальное уравнение Риккати.

Итак, мы получили уравнение продолжения ковариационной матрицы ошибки. Путем интегрирования уравнения (2.3) с начальным условием  получим прогнозируемую оценку ковариационной матрицы ошибки .

2.2 Матрица коэффициентов обратной связи

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

.                                                                               (2.4)

Оценка  базируется на наблюдениях (или измерениях) . Таким образом, она является функцией априори оценки  и измерении  и может быть записана в следующем виде:

.                                                                (2.5)

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

,                                              (2.6)

.

Подставляя в уравнение (2.6) значения  и  из выражений (2.1) и (2.5) и учитывая тот факт, что случайные последовательности  и  некоррелированы, после некоторых преобразований [7], выражаем матрицы коэффициентов обратной связи в следующем виде:

,                                                                               (2.7)

.                                             (2.8)

Эта матрица коэффициентов является функцией от априори значения ковариационной матрицы ошибки.

2.3 Коррекция значения ковариационной матрицы ошибки оценки вектора состояния

Теперь выведем формулу для апостериори значения ковариационной матрицы ошибки (для корректировки ковариационной матрицы ошибки). По определению она записывается в виде

,                                                       (2.9)

где

, .                                   (2.10)

Подставляя (2.7) в (2.5), получим

.                                                 (2.11)

Вычтем  из обеих частей (2.11) и подставим в него значение  в соответствии с уравнением (2.4). В результате получим

или с учетом (2.10)

.                                           (2.12)

Подставляя (2.12) в (2.9) и принимая во внимание, что  получим

.                                                               (2.13)

2.4 Фильтр Калмана для нелинейных систем (расширенный фильтр)

Большинство динамических систем и сенсоров не являются линейными. Но методы фильтрации, разработанные для линейных систем, могут быть применены и для нелинейных систем. В этом разделе рассматриваются методы, которые позволят перенести методологию фильтрации Калмана на нелинейные задачи.

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

Табл.2.2

Модель

Непрерывное время

Дискретное время

Система

Измерения

Применяемый метод линеаризации требует, чтобы функции f и h были дважды непрерывно дифференцируемые. Обозначим символом  малое отклонение от оцениваемой траектории:

,

тогда  в окрестности  может быть представлена в виде

.

Таким образом, получаем

,

где

.                                                       (2.14)

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

 и .

Слагаемые выше первого порядка малости здесь опущены. Если в разложении мы пренебрегаем членами высокого порядка, то возмущение  может быть представлено в виде

,

где

.                                                                     (2.15)

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

,                                                                                                                   (2.16)

.                                                                   (2.17)

2.5 Использование фильтра Калмана

Подведем некоторый итог всего вышеизложенного на примере непрерывно-дискретной системы (непрерывная динамическая система с дискретными измерениями), которая и описывает большинство реальных задач. Такая система записывается следующим образом:

,

.

Допустим, что в момент времени tk-1(+) получены апостериори значения  и , то есть на шаге tk-1(+) задача фильтрации выполнена и теперь необходимо определить  и .

Априори значения оценок вектора состояния и ковариационной матрицы ошибки  и  можно получить путем интегрирования модельного уравнения и уравнения типа Риккати

,

с начальными условиями  и . Так как рассматриваем нелинейную задачу, то для определения матриц  воспользуемся уравнением (2.16). Выполнив эту процедуру, получим  и .

Далее, находим линеаризованную матрицу чувствительности  согласно уравнению (2.15). Определяем матрицу коэффициентов обратной связи , используя уравнение (2.8). Находим апостериори значения  и  по формулам (2.11) и (2.13).

2.6 Проблема адаптивного фильтра. Метод “Уточнения ковариации”

Одной из основных проблем, с которой сталкиваются при работе с фильтром Калмана, является то, что как модель самой системы, так и модель измерений неидеальны. Статистические свойства шумов, а именно, их ковариационные матрицы являются теми параметрами, от которых зависит точность оценивания вектора состояния системы, т. е. точность работы фильтра.

Обычно, начальные значения ковариационных матриц шума измерений и шума процесса выбирается путем анализа некоторых эмпирических данных или путём моделирования различных ситуаций и далее считается постоянной. Если в процессе использования этих данных выясняется, что работа фильтра неудовлетворительна, то требуется новая настройка ковариационных матриц шумов. Появляется задача адаптивной фильтрации – параллельно с оценкой вектора состояния системы оценивать статистические свойства шумов процесса и измерений (ковариационные матрицы) с целью повышения точности работы фильтра. Один из методов такой оценки – метод “уточнения ковариации”. Этот подход базируется на идее, что ковариационная матрица шума процесса равна ковариационной матрице ошибки измерения, то есть . Введём величину, называемую прогнозируемая разность ошибки измерения . Для придания алгоритму большей статистической значимости будем рассматривать выборку из N, сечение случайного процесса .

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

Общую схему для оценки ковариационной матрицы шума процесса R и шума измерений Q можно представить в виде табл.2.3.

Табл.2.3

 

Действия

Оценка R

Оценка Q

1

Начальные условия

j=1,

j=1,

2

Аппроксимация шума

3

Вычисление

4

j=N

Подставляем  в фильтр

Подставляем  в фильтр

5

j<N

Увеличиваем j на 1 и переходим к шагу 2

3. Уравнения движения макета

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

3.1. Анализ сил и моментов, действующих на макет

На макет действуют в основном следующие силы и моменты:

·       сопротивление движению со стороны воздуха;

·       реакция шарнира (опоры), которая состоит из момента реакции и силы реакции, приложенной к точке опоры;

·       сила тяжести.

Сопротивление воздуха не поддается простому описанию, тем более что оно зависит от аэродинамических характеристик самого макета. Из-за малых линейных скоростей точек макета пренебрегаем им.

Оценим минимальный диаметр d струны подвеса. Механическое напряжение не должно превышать напряжение σт текучести материала подвеса с коэффициентом запаса n, то есть

Здесь M – масса макета. Примем M=1 кг и коэффициент запаса 1.5. В случае стали (σт= 200 МПа) диаметр составит d>0.35 мм. В случае сплава титана (σт=900 МПа) d>0.17 мм.

Определим момент реакции шарнира в случае подвеса как упругий момент закрученной нити. Погонный угол закрутки φ=2πk/L, где k – число оборотов нити, L – ее длина. Механический момент составит [8]

                                                                           (3.1)

где G=8·1010 Па - модуль Юнга второго рода для стали, Ip=πD4/32 - полярный момент инерции круглого стрежня, f – модуль кручения.

3.2. Вывод уравнений движения тела на струне

Рис.3.1. Системы координат

Рис.3.2. Силы, действующие на тело

Сделаем предположение, что струна бесконечно длинная. Как следствие этого считаем, что точка подвеса всегда находится в одной и той же горизонтальной плоскости. Таким образом, положение точки подвеса A можно задать двумя переменными x и y (рис.3.1) в неподвижной системе координат Oxyz [9]. Еще три переменные (углы Эйлера φ, ψ и θ) задают ориентацию макета. Таким образом, получаем пятимерный вектор состояния макета. Свяжем с макетом систему координат Ax`y`z` как показано на рис.3.1. Пусть центр масс макета лежит на оси z` и смещен вниз на расстояние d от точки подвеса. Его координаты в системе координат Ax`y`z` суть (0,0,-d). К центру масс C тела приложена сила тяжести Fт (рис.3.2). Из сделанного предположения следует также, что струна всегда вертикальна. Следовательно, сила реакции N подвеса направлена всегда вверх. Так как точка A подвеса всегда лежит в одной и той же горизонтальной плоскости, то сила реакции N не совершает работы.

Введем обозначения: A, B и C суть его главные моменты инерции. Пусть, для простоты, эти главные оси инерции оси направлены параллельно осям координат x`, y` и z`.

Составим уравнения Лагранжа второго рода. Кинетическая энергия тела складывается из энергии поступательного движения центра масс и энергии вращения вокруг центра масс. Найдем скорость движения центра масс. Она складывается из переносной скорости движения точки подвеса и относительной скорости движения центра масс вокруг точки подвеса [10]. В проекциях на оси x, y, z:

Кинетическая энергия, согласно теореме Кенига, запишется в виде

где p, q, r - проекции угловой скорости вращения макета на его главные оси инерции, вычисляемые по формулам

Потенциальная энергия складывается из потенциальной энергии макета в поле силы тяжести и потенциальной энергии закрученной струны. Таким образом, можно записать

.

Лагранжиан системы L=T-П. Уравнения Лагранжа второго рода запишутся в координатном виде

(3.2)

Очевидно, что vCx и vCy являются первыми интегралами системы, они выражают закон сохранения импульса.

Теперь будем считать углы θ и φ малыми, а центр масс макета будет находиться в покое. Тогда, как видно из уравнений движения (3.2), получаем одномерное движение

.                                                                                     (3.3)

Это – уравнение крутильных колебаний, которое будет использовано для задания модели движения в фильтре Калмана.

4. Применение фильтра Калмана

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

4.1 Перманентное вращение

В первой модели предполагается, что тело на струне вращается с постоянной скоростью. Понятно, что на самом деле тело меняет скорость своего вращения вследствие силы упругости струны, но, как показывает практика, если струна достаточно длинная, такая модель неплохо описывает движение системы на небольших промежутках времени. Кроме того, фильтр, основанный на этой модели, демонстрирует насколько “страдает” оценка состояния системы, если выбрана достаточно грубая динамическая модель. Итак, пусть движение описывается уравнением . Вектор состояния и динамическое уравнение в нашем случае

,      ,

Уравнение можно записать по-другому

.                                                                       (4.1)

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

.

Проинтегрируем уравнение (4.1), получим

         .

Пусть ковариационная матрица ошибки имеет следующие элементы:

.

Тогда матричное уравнение Риккати (3) с начальными условиями  имеет следующее решение:

            .            (4.2)

Матрица коэффициентов обратной связи, согласно (2.8), имеет вид

.                                                                                    (4.3)

Тогда апостериори оценка вектора состояния и матрицы ковариации из (2.11) и (2.13) соответственно, имеют вид

,                                                   (4.4)

.                                                     (4.5)

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

4.2 Крутильные колебания

Используем полученное уравнение движения макета (3.3), которое запишем так: . Здесь I – момент инерции тела относительно оси вращения,  - угол поворота тела относительно положения равновесия. Обозначим . Вектор состояния и динамическое уравнение имеют вид

,    

или

.                                                                  (4.6)

Модель измерения имеет вид

.

Проинтегрируем уравнение (4.6) и получим его решение

, .

Матричное уравнение Риккати (2.3) с начальными условиями  имеет решение

,                                 (4.7)

.

Коэффициент обратной связи согласно (2.8) и оценка вектора состояния и матрицы ковариации из (2.11) и (2.13) будут иметь такой же вид, как и для модели перманентного движения: (4.3), (4.4) и (4.5) соответственно

5. Математическое моделирование

определения ориентации

Проведём математическое моделирование определения ориентации тела, подвешенного на струне, имитируя зашумлённые показания датчика угловой скорости, основанные на более реальной модели движения – крутильные колебания с затуханием. Программа математического моделирования реализована в интегрированной среде МATLAB.

5.1 Структура программы

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



На вход программы подаются с некоторой частотой имитируемые измерения угловой скорости (блок-схема программы приведена на рис.5.2). Эти данные усредняются с заданным параметром усреднения в целях уменьшения ошибки поступающих данных. Усреднённая угловая скорость поступает на вход сразу трём программным модулям: интегрирование кинематического уравнения, Фильтр Калмана (ФК) модель перманентного вращения, ФК модель крутильных колебаний. На выходе соответственно имеем: результат интегрирования, апостериори оценки угла отклонения макета от положения равновесия и угловой скорости двух ФК. Все эти значения и значение усреднённой угловой скорости подаются на вход модуля программы “Уточнение ковариации”, который работает по алгоритму, описанному в табл.2.3. Он с некоторым периодом подставляет в фильтры Калмана оцененное значение матрицы шума системы Q. Все выходные значения программы записываются в файл. В конце каждого блока выходной информации записывается время с начала работы программы.

После отработки программы строятся графики на основе данных из файла, содержащего результаты вычислений (рис.5.3) Первый график отображает результат интегрирования кинематического уравнения и оценки двух фильтров Калмана угла отклонения от положения равновесия. На втором графике отображены показания датчика и оценки фильтров относительно угловой скорости. Третий и четвёртый графики отображают покомпонентную разницу между вектором фазового состояния и выходом двух фильтров. Эти графики необходимы, что иметь представление о точности оценок фильтров Калмана.

5.2 Результаты моделирования

Сначала моделировалась работа неадаптивного фильтра Калмана, который работает по такой же схеме, как на рис.5.2, но без звена уточнения матрицы ковариации шума системы. Результаты представлены на рис.5.3.

Как видно из графиков, фильтры Калмана, основанные на неточных моделях процесса, начинают подстраиваться, но подстраиваются весьма плохо, поэтому ошибка в определении вектора состояния достаточно велика. Причину плохой оценки можно понять, если взглянуть на графики элементов матрицы коэффициентов матрицы ковариации (рис.5.4) и на выражение (4.4). Эти графики определяются только решением матричного уравнения Риккати (2.3) для матрицы ковариации, которое имеет вид (4.2) и (4.7). Согласно этому решению элементы матрицы очень быстро сходятся к нулю, и, как можно судить из выражения (4.4), оценка компонент вектора состояния становится слабо чувствительной к измерениям, и строится только на основании заложенной математической модели. Это было бы не плохо, если математическая модель очень точно описывает движение тела, однако для всех реальных задач это не так. Поэтому фильтр Калмана плохо справляется со своей задачей – оценкой вектора состояния системы по имеющимся зашумленным измерениям. Хотя, как можно увидеть, оценка ориентации фильтром, основанным на модели перманентного вращения, достаточно точна. Это обусловлено особенностью решения уравнения Риккати (см. рис.5.4). Благодаря такому решению, фильтр остался чувствительным к измерениям.

Итак, из-за неточности математических моделей не следует использовать фильтр Калмана в таком виде, необходимо отслеживать “пригодность” фильтра. Сделаем это с помощью алгоритма “уточнение ковариации”, описанного в табл.2.3. Зададим размер выборки, на основе которой оценивается матрица ковариации, равным N=1000. Результаты представим на рис.5.5.

Как видно из графиков, точность фильтров значительно улучшилась. Это произошло благодаря тому, что раз в 10 секунд в фильтры Калмана подставлялись оцененные элементы матрицы ковариации (см. рис. 5.6), чувствительность к измерениям резко возрастала, что следует из выражения (4.4), и оценка вектора состояния становилась приемлемой, подстраиваясь под реальные движения тела. Следующие десять секунд элементы матрицы ошибки вычисляются по решению уравнения Риккати с подставленными начальными условиями, и, как видно из рис. 5.6, снова быстро сходятся к нулю. Оценки фильтров снова начинают зависеть практически только от вложенной в них модели и расходиться с реальным движением.

Возникает закономерный вопрос: почему бы не сделать элементы матрицы ковариации постоянными? Тогда фильтры будут всегда давать оценки, максимально приближенные к реальным движениям.

Дело в том, что в этом случае, оценки фильтров Калмана мало будут отличаться от шума датчика, никакой плавной кривой не получится, задача фильтра – устранение шума – не будет выполняться. Поэтому тут нужно найти “золотую середину” между целью иметь максимально близкие к реальности оценки и целью устранить шум датчика. Это можно сделать, варьируя размер выборки N. Случай N=1 соответствует крайнему варианту – фильтр не устраняет шум. При N>2000 и частоте входных параметров 100 Гц оценки фильтра становятся неудовлетворительными. Для рассматриваемой задачи определения ориентации тела, подвешенного на струне, оптимум составляет N=100 (см. рис.5.7).


Графики ошибки фильтров, особенно для ошибки скорости, являются только лишь оценочными (см. нижний правый прямоугольник на рис.5.7,). Здесь оценки фильтров сравниваются с зашумлёнными измерениями датчика. То, что ошибка измерения колеблется около нуля, говорит о “хорошей” работе фильтров. Приемлемость оценок фильтра можно условно определить следующим образом: оценки фильтра должны находиться в области шума. На рис.5.8 это условие полностью достигается фильтром, основанным на модели крутильных колебаний, и частично фильтром перманентного вращения. Поэтому можно сказать, что первый фильтр даёт более точную оценку.


6. Проведение измерений с использованием макета

Рассмотрим схему макета, собранного в ИПМ им.М.В.Келдыша РАН, устройство измерительных элементов, установленных на макете, и результаты проведенных испытаний.

6.1 Схема макета

Макет представляет собой автономный прибор, включающий в себя систему электропитания, систему определения ориентации и систему связи с стационарным компьютером. Схема макета изображена на рис.6.1, сам макет изображен на рис.6.2.

Система электропитания включает в себя два аккумулятора и преобразователь напряжения. Преобразователь формирует из 17 вольт от аккумуляторов 5 вольт, необходимых для питания бортового компьютера и микромеханического гироскопа iMEMS ADIS16100.

Рис.6.1. Схема макета

Система определения ориентации состоит из бортового компьютера Jrex-CE (далее БК), оптоволоконного гироскопа VG-910D и микромеханического гироскопа iMEMS ADIS16100. БK производит считывание информации об угловой скорости с одного из гироскопов, обрабатывает её с помощью программного обеспечения (далее ПО) в режиме реального времени и сохраняет результаты в файл. Передача информации производится через USB порты БК, используется разветвитель USB-Hub.

Система связи состоит из устройства Bluetooth и БК. Устройство Bluetooth позволяет обмениваться файлами со стационарным компьютером, к которому также подключено устройство Bluetooth. Посредством обмена файлами происходит считывание с БК результатов вычислений, а также управление ПО БК.

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

6.2 Характеристика и устройство датчиков угловой скорости

Волоконно-оптический гироскоп VG-910D

VG-910D – это одноосный прецизионный датчик вращения, выполненный по волоконно-оптической технологии, производства компании ФИЗОПТИКА [11]. Его функционирование основано на релятивистском эффекте Саньяка. Датчик содержит волоконно-оптический преобразователь и блок электроники, реализующий современные алгоритмы обработки аналоговых сигналов. Он предназначен для использования в различных областях техники с целью измерения угловой скорости вращения. Выходной цифровой сигнал пропорционален угловой скорости с некоторым масштабным коэффициентом. Диапазон измерений град/сек, цифровой выход угловой скорости – 24-х разрядное число в дополнительном коде, входное напряжение – 5V, может питаться через USB кабель. Датчик весит 130 грамм при габаритах 80мм19мм. Частота выходных данных – 300Гц.

VG-910D состоит из двух основных узлов (рис.6.3):

·   волоконно-оптический преобразователь – кольцевой интерферометр, содержащий волоконный контур (катушка), сварные волоконно-оптические ответвители, волоконно-кристаллический поляризатор, пьезокерамический модулятор и светодиодный излучательный модуль;

·   плата обработки с аналого-цифровым преобразователем.

Рис.6.3. Схема датчика VG-910D

Принцип действия датчика основан на зависимости времени распространения света по замкнутому вращающемуся контуру от направления обхода. В интерферометре излучение делится на два луча равной интенсивности, которые распространяются в контуре в противоположных направлениях, а затем смешиваются и интерферируют. При вращении между волнами возникает разность фаз (фаза Саньяка), которая и является мерой угловой скорости. Коэффициент пропорциональности между угловой скоростью и фазой Саньяка называется оптическим масштабным коэффициентом преобразователя. Его величина определяет диапазон измеряемых угловых скоростей.

Интегральный гироскоп на базе iMEMs-технологии

ADIS16100 – интегральный гироскоп на базе технологии iMEMS (integrated Micro Electro Mechanical System), устройство фирмы Analog Devices, интегрирующее на одном кремниевом кристалле датчик угловой скорости и электронику, обеспечивающую формирование и предварительную обработку сигнала [12]. Гироскоп изображен на рис.6.4.

Несмотря на меньшую в сравнении с прочими гироскопами точность, микромеханические гироскопы iMEMS обладают целым рядом уникальных достоинств. Прежде всего – это малые габариты и масса. Датчик угловой скорости ADIS16100 выпускается в миниатюрном корпусе размером 773мм, вес такого прибора не превышает 0,5 г. Диапазон измерения датчика град/c, цифровой выход угловой скорости – 12-ти разрядное число. Частота выходных данных – 1400Гц.

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

6.3 Программное обеспечение макета

Программное обеспечение (ПО) – это программа, работающая на БК макета, осуществляющая считывание информации с датчика угловой скорости, её обработку и сохранение результатов вычислений в файл. Программа написана в среде C++ Builder 6.

6.3.1 Работа программы

С датчика поступает цифровой сигнал, который преобразуется калибровкой в угловую скорость (рис.6.5). Значения угловой скорости за несколько измерений суммируются, после чего эта сумма делится на количество измерений, таким образом, на выходе получаем среднеарифметическое значение угловой скорости. Усреднение производится для того, чтобы уменьшить шум измерений и уменьшить число вычислений программы в единицу времени. У ПО есть верхний предел по частоте выходных данных (100Гц). Выше этого значения программа не успевает обрабатывать поступающие с датчика данные и обработка в режиме реального времени становится невозможной. Поэтому при установке параметров работы программы нужно отдельно обратить внимание на усреднение: зная частоту подачи данных устанавливать такое усреднение (то есть выбирать количество измерений для усреднения), чтобы частота выходных данных программы не превышала предела.

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

Все возможные настройки программы при её запуске загружаются из файла “Parameters.txt”, который находится в той же папке, что и программа. Во время эксперимента нет доступа к интерфейсу программы, работающей на БК, поэтому управление программой производится с помощью Bluetooth посредством изменения содержимого файла “Parameters.txt.

6.3.2 Интерфейс программы

Главное окно программы содержит две области (рис.6.6):

·       Settings – установка рабочего датчика, выбор порта, через который будет поступать информация, и установка усреднения. Задав эти параметры программы можно нажать кнопку “Open port”, и автоматически пойдет процесс вычисления

·       Some information – содержит информацию о частоте входных значений с датчика, о выходной частоте программы, отображает время работы программы и число принятых байт в одной посылке.

Главное окно содержит меню со следующими разделами:

·       File – выход из программы

·       Sensor calibration – установка калибровочных коэффициентов для датчиков VG-910D и ADIS16100.

·       Initial conditions – задание начального вектора состояния макета.

·       Model of measuring – задание параметров модели измерений (квадрата ошибки измерений  и значений элементов матрицы ошибки для ФК).

·       Graphics – вывод в режиме реального времени графиков угловой скорости и угла отклонения от положения равновесия. Графики отображают все выходные данные программы.

·       Output values – вывод на экран выходных значений двух фильтров и результата интегрирования.

6.4 Результаты экспериментов

Были проведены эксперименты с использованием макета тела, подвешенного на струне. Результат одного из экспериментов приведён на рис.6.7.

Как видно из графиков, реальное движение макета – это крутильные колебания с довольно сильным затуханием. Кроме того, под воздействием каких-то внешних факторов (например, малое колебание точки подвеса макета), подвешенное тело начинает двигаться по не совсем понятному закону (рис.6.8). Такие импульсные и случайные возмущения трудно учесть и ввести в математическую модель движения.

Однако фильтры Калмана и в этом случае, когда движение вообще не соответствует введенной в них математической модели, даёт неплохую оценку состояния тела.

Если рассмотреть ближе импульсное возмущение (рис.6.9), то видно, что с помощью оценки фильтров мы имеем достаточно точные представления о фазовом состоянии макета в короткий промежуток времени действия возмущения.

Зачастую, во время проведения эксперимента не удавалось полностью обеспечить одномерное движение макета, оставались малые колебания относительно подвижной горизонтальной оси, которые тоже давали возмущения в движение тела. Но, как видно из рис.6.10, и в этом случае фильтры Калмана давали достаточно точное представление о реальном движении макета.

Во время экспериментов выяснилось, что из-за физических свойств нити нельзя достаточно точно определить положение равновесия макета – конкретной точки, находясь в которой, макет будет находиться в состоянии покоя. Получалась целая довольно обширная область таких точек. Таким образом, определялась ориентация макета относительно наперед неопределённого положения. Проблема состоит в том, что в этом случае очень неточно задаются начальные условия для математических моделей, на которых построены фильтры, а это особенно важно для модели вращательных колебаний, так как действующий на макет момент со стороны упругой нити пропорционален углу отклонения от положения равновесия. Ожидалось, что фильтр будет давать неудовлетворительную оценку. Но, как показал эксперимент, работа фильтра не сильно пострадала из-за этого, казалось, принципиального момента (рис. 6.7). Даже, когда “промах” в определении начального положения составлял около ста градусов, оценка вектора состояния была достаточно точной.

7. Заключение

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

·                расширенный адаптивный фильтр Калмана, основываясь на грубой математической модели, достаточно точно оценивает вектор состояния тела,

·                фильтр Калмана даёт точные оценки и в случае возмущающих воздействий, действующих на тело,

·                оценка также не страдает, если неточно заданы начальные условия, от которых сильно зависит математическая модель фильтра.

Таким образом, фильтр Калмана, основываясь на зашумлённых измерениях датчика угловой скорости и неточной математической модели движения тела, даёт достаточно точные оценки вектора фазового состояния тела, которые можно использовать, к примеру, для управления ориентацией. Однако точность рассмотренного способа определения ориентации можно улучшить, если использовать в паре с показаниями датчика угловой скорости измерения датчика направления, например, солнечного датчика.

8. Благодарности

Работа выполнена при частичной поддержке РФФИ (грант 06-01-00389) и DAAD (программа Leonhard Euler, реферат 325).


9. Литература

1.                   А.М. Овчинников, М.Ю. Овчинников, А.С. Середницкий. Лабораторный стенд для отработки алгоритмов определения движения по снимкам звездного неба/ Препринт ИПМ им. М.В. Келдыша РАН. – М., 2006. – №43.

2.                   А.А. Дегтярёв, М. Грасси, А. Перрота, М.Ю. Овчинников. Тестирование и калибровка миниатюрного солнечного датчика, разработанного на основе APS-технологии / Препринт ИПМ им. М.В. Келдыша РАН. – М., 2005. – №94.

3.                   А.А. Ильин , Н.В. Куприянова, М.Ю. Овчинников , В.И. Пеньков, А.С. Селиванов. Анализ вращательного движения первого российского наноспутника ТНС-0 по результатам летных испытаний/ Препринт ИПМ им. М.В. Келдыша РАН. – М., 2006. – №18.

4.                   М Л. Пивоваров. Определение ориентации ИСЗ с использованием измерений угловых скоростей // Космические исследования. – 1985. – Т. 23, вып. 3. – С. 331-334.

5.                   В.Н. Бранец. О точности решения кинематических уравнений // Космические исследования. 1982. – Т. 20, вып. 2. – С. 184-190.

6.                   С.В. Барабаш, И. Ю. Кирюшкин, У. Норберг, М.Ю. Овчинников. Методы управления углового движения нано-спутника Munin/ Препринт ИПМ им. М.В. Келдыша РАН. – М., 2005. – №94.

7.                   А.А. Дегтярёв, Ш. Тайль. Элементы теории адаптивного расширенного фильтра Калмана / Препринт ИПМ им. М.В. Келдыша РАН. – М., 2003. – №26. – 35 с.

8.                   А.Д. Гладун. Лабораторный практикум по общей физике/М.: МФТИ, 2004. – 316 с.

9.                   М.Ю. Овчинников, Е.А. Цветков. Проектирование имитатора геомагнитного поля в составе лабораторного стенда для отработки способов управления ориентацией микроспутников / Препринт ИПМ им. М.В. Келдыша РАН. – М., 2005. – №55. – 29 с.

10.              А.Ю. Ишлинский, В.А. Стороженко, М.Е. Темченко. Вращение твёрдого тела на струне и смежные задачи / М.: Наука, 1991. – 330 с.

11.              http://www.fizoptika.ru

12.              http://www.mems.ru/