CC BY
246
МЕХАНИКА
УДК 531 (075.8); 629.7.05(075)
ИССЛЕДОВАНИЕ АЛГОРИТМОВ ОПРЕДЕЛЕНИЯ ИНЕРЦИАЛЬНОЙ ОРИЕНТАЦИИ ДВИЖУЩЕГОСЯ ОБЪЕКТА
Ю. Н. Челноков1, С. Е. Переляев2, Л. А. Челнокова3
1 Челноков Юрий Николаевич, доктор физико-математических наук, профессор кафедры математического и компьютерного моделирования, Саратовский национальный исследовательский государственный университет имени Н. Г. Чернышевского; заведующий лабораторией механики, навигации и управления движением Института проблем точной механики и управления РАН, Саратов, ChelnokovYuN@info.sgu.ru
2 Переляев Сергей Егорович, доктор технических наук, главный научный сотрудник, ООО «Аэроспецпроект», Жуковский, Московская обл., Pers2030@yandex.ru
3 Челнокова Людмила Александровна, младший научный сотрудник лаборатории механики, навигации и управления движением Института проблем точной механики и управления РАН, Саратов, ChelnokovaLA@mail.ru
Исследуются новые и известные алгоритмы работы бесплатформенных инер-циальных навигационных систем, предназначенные для высокоточного определения параметров ориентации движущегося объекта (параметров Родрига-Гамильтона (Эйлера)) в инерциальной системе координат. В основе построения новых алгоритмов лежит использование классического кватерниона поворота Гамильтона, кватерниона с нулевой скалярной частью, который ставится в соответствие классическому кватерниону поворота с помощью кватернионного аналога формулы Кэли, а также используется новое кватернионное дифференциальное уравнение инерциальной ориентации движущегося объекта. Новые алгоритмы построены с помощью метода последовательного приближения Пикара. В этих алгоритмах в качестве входной информации используется интегральная первичная информация об абсолютном угловом движении объекта. Показывается, что новые алгоритмы имеют преимущества перед известными алгоритмами аналогичного порядка в смысле точности и сложности.
Ключевые слова: движущийся объект, параметры Родрига-Гамильтона (Эйлера), инерциальная ориентация, кватернион Гамильтона, кватернионная матрица, формула Кэли, четырехмерный кососимметрический оператор.
DOI: 10.18500/1816-9791 -2016-16-1-80-95 ВВЕДЕНИЕ
В настоящей работе для нахождения параметров инерциальной ориентации движущегося объекта (твердого тела) с помощью бесплатформенной инерциальной навигационной системы (БИНС) используется дифференциальное нелинейное кватернионное или матричное кинематическое уравнение типа Риккати [1]. Переменная в этом уравнении — кватернион с нулевой скалярной частью (ассоциированный кватернион) или четырехмерная кососимметрическая
матрица, которые ставятся в соответствие классическому кватерниону конечного поворота или ква-тернионной матрице поворота с помощью формул Кэли. Эта переменная отвечает четырехмерному кососимметрическому оператору, который соответствует трехмерному вектору конечного поворота в = tg(</4)e движущегося объекта (< — эйлеров угол поворота, е — единичный вектор эйлеровой оси поворота объекта).
В статье рассмотрены новые кватернионный одношаговый алгоритм третьего порядка точности и кватернионные двухшаговые алгоритмы третьего и четвертого порядков точности в четырехмерных кососимметрических операторах, предназначенные для вычисления параметров инерциальной ориентации объекта. Алгоритмы построены с помощью метода последовательного приближения Пикара, используют в качестве входной информации интегральную первичную информацию измерителей абсолютной угловой скорости объекта (приращения интегралов от проекций вектора абсолютной угловой скорости объекта на связанные с ним координатные оси (квазикоординаты)). С помощью математического моделирования показывается, что эти алгоритмы имеют преимущества перед известными алгоритмами аналогичного порядка в смысле точности и сложности. Статья является развитием работ [2-5].
1. КВАТЕРНИОННЫЕ КИНЕМАТИЧЕСКИЕ УРАВНЕНИЯ ОРИЕНТАЦИИ ТИПА РИККАТИ
Ортогональная кватернионная матрица поворота п = п(А) [1], сопоставляемая кватерниону поворота Гамильтона А, связана с четырехмерной кососимметрической матрицей к формулами Кэли [1,6]:
п = (Е - к)(Е + к)
-1
п(А) =
к= (Е — п)(Е + п)- 1,
/Ао —А1 —А2 —Аз \ /0 —к1 —к2 —к3\
А1 Ао А3 —А2 , к = к1 0 к3 —к2
А2 —А3 Ао А1 к2 —к3 0 к1
А3 А2 —А1 Ао к3 к2 —к1 0
(1) (2)
(3)
п = п(А) = п(Ао + А„) = п(сов(</2) + е э1п(</2)), е5 = ех = е^ + в212 + ез1з.
Здесь Е — четырехмерная единичная матрица, запись вида п(А) означает, что кватернионная матрица п сопоставляется кватерниону А, записанному в круглых скобках; ех и е^ — отображения единичного вектора е на связанный с объектом базис X и инерциальный базис £ соответственно, е^ — проекции вектора е на оси систем координат X и £ (одинаковые при условии, что в начальном положении одноименные оси систем координат X и £ совпадали); ^, 12, — векторные мнимые единицы Гамильтона.
Кососимметрическая матрица к имеет вид
к = ^(</4)
0 —е1 —е2 —е3
е1 0 е3 —е2
е2 —е3 0 е1
е3 е2 —е1 0
(4)
Отметим, что если трехмерной кососимметрической матрице, связанной с матрицей направляющих косинусов формулой Кэли, соответствует вектор к = —tg(</2)e, то четырехмерной кососимметриче-ской матрице к, связанной этой формулой с кватернионной матрицей поворота п, соответствует, как это следует из (4), вектор к = —tg(</4)e. При этом если трехмерная кососимметрическая матрица и соответствующий ей вектор к не определены для углов < = ±п, то четырехмерная кососимметрическая матрица к и соответствующий ей вектор к не определены для углов < = ±2п.
Матричное кинематическое уравнение, использующее в качестве переменной четырехмерную ко-сосимметрическую матрицу к, определяемую формулой (4), может быть записано в одной из двух следующих форм, имеющих вид матричных уравнений Риккати [1]:
4Г = —(Е + к)0(Е — к),
(5)
4к* = кОк + Ок - кО - О.
(6)
В этих уравнениях верхняя точка — символ дифференцирования по времени, О — четырехмерная кососимметрическая матрица, сопоставляемая отображению шх вектора ш абсолютной угловой скорости объекта на связанный базис X, имеющая вид
О =
(0 —ш1 —Ш2 —Шз\
Ш1 0 Шз —Ш2
Ш2 —Шз 0 Ш1
\шз Ш2 —Ш1 0
(7)
где ш, — проекции вектора ш на оси связанной системы координат X.
Отметим, что в задачах определения ориентации движущегося объекта используется также [7,8] матричное кинематическое уравнение, где в качестве переменной берется трехмерная кососиммет-рическая матрица к, получаемая вычеркиванием первой строки и первого столбца в четырехмерной матрице к, определяемой второй из формул (3). Это уравнение формально отличается от уравнения (6) числовым коэффициентом, стоящим перед производной от используемой переменной (вместо коэффициента 4 стоит коэффициент 2):
2к* = кОк + Ок - кО - О.
Матрица О в этом случае — трехмерная кососимметрическая матрица угловых скоростей, получаемая вычеркиванием первой строки и первого столбца в четырехмерной матрице О, определяемой формулой (7), а трехмерная кососимметрическая матрица к связана с матрицей с направляющих косинусов углов формулой Кэли [1]: с = (Е — к)(Е + к)-1.
Кососимметрические матрицы третьего (нечетного) и четвертого (четного) порядков имеют качественно различные свойства: если кососимметрические матрицы третьего порядка являются особыми (их определители равны нулю), то кососимметрические матрицы четвертого порядка — нет (их определители всегда отличны от нуля). Кроме того, если многочлен любой степени от кососимметрической матрицы третьего порядка сводится к многочлену второй степени, то многочлен любой степени от кососимметрической матрицы четвертого порядка сводится к многочлену первой степени. Последнее обстоятельство делает использование матричного кинематического уравнения (6) в четырехмерных кососимметрических операторах при построении алгоритмов определения ориентации движущихся объектов с помощью БИНС более эффективным по сравнению с использованием кинематического уравнения в трехмерных кососимметрических операторах.
В кватернионной записи кватернионные аналоги формул Кэли, связывающие кватернион поворота Л с соответствующим трехмерным вектором к (точнее, с кватернионом к, имеющим нулевую скалярную часть), имеют вид
Л = (1 + к)-1 о (1 — к) = (1 — к) о (1 + к)-1, (8)
к = (1 + Л)-1 о (1 — Л) = (1 — Л) о (1 + Л)-1, (9)
где
Л = Ао + А111 + Л2Ь + Аз ¡з = Ао + Л„ = соб(^/2) + вт(р/2)(е111 + в212 + еэ1э), к = к111 + к212 + кз1з = — tg (^/4)(в111 + в212 + ез 1з),
центральный кружок — символ кватернионного произведения.
Из этих формул следует, что четырехмерной кососимметрической матрице к, порождаемой формулой Кэли, отвечает кватернион к = —tg (^/4)(е111 + е212 + ез 1з) с нулевой скалярной частью.
Кватернионные кинематические уравнения типа Риккати, соответствующие матричным уравнениям (5) и (6), имеют следующий вид [1]:
4к* = —(1 — к) о шх о (1 + к), (10)
4к = к о шх о к + к о шх — шх о к — шх = к о шх о к — 2шх х к — шх, (11)
где кватернионная переменная к = к111 + к2Ь + к313 и кватернион угловой скорости шх = о11 + + о2Ь + о3Ь — кватернионы с нулевыми скалярными частями, х — символ векторного произведения.
Отметим, что векторы конечных поворотов в = tg(</4)e и в = С^(</4)е и соответствующие им векторные кинематические уравнения рассматривались в [9].
Кроме того, отметим аддитивное вхождение матрицы угловых скоростей О в матричное кинематическое уравнение (6), кватерниона угловой скорости шх — в кватернионное кинематическое уравнение^), а вектора угловой скорости ш — в часто используемое векторное кинематическое уравнение для вектора ^ = <е эйлерова конечного поворота твердого тела. Поэтому первое приближение решений этих уравнений на шаге интегрирования, построенное методом последовательного приближения Пикара, есть приращение интеграла на шаге интегрирования от вектора абсолютной угловой скорости объекта, непосредственно измеряемое на борту движущегося объекта большинством современных датчиков угловой скорости. Это делает указанные уравнения более удобными для построения алгоритмов определения инерциальной ориентации движущегося объекта, использующих интегральную первичную информацию об угловом движении объекта, по сравнению с классическим кватернионным кинематическим уравнением:
2А^ = А о шх,
использующим в качестве переменной кватернион ориентации А.
Известно, что компонента АО скалярной части и компоненты А* (г = 1, 2,3) векторной части кватерниона приращения А* кватернионной переменной А, вычисляемые на шаге интегрирования, отличаются на несколько порядков: величина А* близка к единице, в то время как величины А* имеют порядок, равный 10-5 и меньше (в зависимости от вида углового движения объекта и шага интегрирования). Близость компоненты А* приращения А* кватернионной переменной А, вычисляемого на шаге интегрирования, к единице может приводить к потере точности определения ориентации при малой ограниченной разрядной сетке бортового вычислителя. Поэтому при построении алгоритмов численного интегрирования классического кватернионного кинематического уравнения целесообразно вместо переменной АО использовать новую переменную, равную 1 — А*. Во многих алгоритмах, предложенных зарубежными и отечественными авторами, величины А* вычисляются через промежуточные кинематические параметры <, являющиеся проекциями вектора ^ = <е эйлерова конечного поворота объекта на связанные координатные оси. Нами предлагается вычислять величины А* через другие промежуточные кинематические параметры к, являющиеся проекциями вектора к = —tg(</4)e конечного поворота объекта на связанные координатные оси.
2. АЛГОРИТМЫ ОПРЕДЕЛЕНИЯ ИНЕРЦИАЛЬНОЙ ОРИЕНТАЦИИ
В основе решения задач автономного определения ориентации движущегося объекта лежит численное интегрирование дифференциальных уравнений инерциальной ориентации объекта на его борту. Для интегрирования этих уравнений могут быть использованы стандартные методы численного интегрирования дифференциальных уравнений, такие как метод Рунге-Кутта, метод прогноза и коррекции и другие стандартные методы, использующие мгновенную информацию о проекциях о вектора ш абсолютной угловой скорости объекта на связанные координатные оси, а также численные методы, специально разработанные для интегрирования уравнений инерциальной ориентации и использующие либо мгновенную информацию о проекциях о, либо интегральную информацию (приращения интегралов от проекций о угловой скорости объекта на связанные координатные оси).
Выбор того или иного метода интегрирования зависит от того, для какого класса движущихся объектов он предназначен (зависит от законов и параметров ожидаемых вращательных движений объекта). Критериями при выборе метода служат требуемая точность вычисления параметров ориентации и объем необходимых вычислений. Ниже приводятся новые алгоритмы [2-5], построенные методом последовательного приближения Пикара с использованием формул Кэли и кватернионного кинематического уравнения типа Риккати (11). Отметим, что порядок точности алгоритма, например второй, означает, что методическая накапливающаяся погрешность алгоритма пропорциональна Л2, где к — шаг интегрирования.
Алгоритмы инерциальной ориентации движущегося объекта (как правило, алгоритмы вычисления кватерниона Л инерциальной ориентации) строятся по общей рекуррентной схеме, имеющей в кватернионной записи следующий вид:
Л(и) = Л(г„_1) о л*, (12)
3 3
Л = Ао + Лу = Ао + ^ Afcifc, Л* = А* + Л* = А* + ^ Aklk • k=1 k=1
Здесь A* — приращения параметров Эйлера (Родрига - Гамильтона) Aj на шаге интегрирования h = tn — tn-1; Л,, Л: — векторные части кватернионов Л и Л*.
Приращения A*,A*(i = 1, 2,3) (или А*, Л*) параметров Эйлера вычисляются на каждом шаге интегрирования по соотношениям, вид которых зависит от выбранного метода численного интегрирования. Часто, как уже отмечалось, приращения А* вычисляются через промежуточные кинематические параметры p^, являющиеся проекциями вектора ^ = pe эйлерова конечного поворота объекта на связанные координатные оси, по формулам
A* = cos(p/2), Ak = p-1pk sin(p/2) (k = 1, 2, 3), p = (p2 + p2 + p2)1/2. (13)
В свою очередь, переменные p^ находятся по одному из выбранных алгоритмов численного интегрирования (с нулевыми начальными условиями) матричного дифференциального уравнения Стьюль-пнагеля [7,8] (J. Stuelpnagel):
= О — ^(Ok — kO) + -1(1 — pctg (p/2))(k2О + Ok2 — 2k0k), 2 p2 2
в котором k и О — трехмерные кососимметрические матрицы, сопоставляемые векторам ^ и ш соответственно (построены из проекций этих векторов на связанные координатные оси), или эквивалентного ему векторного дифференциального уравнения Борца [10] (J. E. Bortz):
= ш + 1 х ш + -1(1 — Pctg (p/2))<p х (<р х ш). (14)
2 p2 2
Каждое из этих уравнений является существенно нелинейным и содержит особую точку p = ±2п.
Отметим, что при построении многих предложенных алгоритмов ориентации приращение вектора ^ на шаге интегрирования вычисляется с помощью интегрирования уравнения Борца (14), в котором отбрасывается или упрощается [11, 12] третье нелинейное слагаемое, содержащее двойное векторное произведение. Так, в [11] в третьем слагаемом уравнения (14) полагается (1/p2)(1 — (p/2)ctg (p/2)) ^ 1/12, а затем получающееся в этом уравнении выражение (1/2)^ х ш + (1/12)^ х (^ х ш) упрощается. Такие упрощения допустимы лишь при построении алгоритмов, не выше третьего порядка точности. Нами приращения параметров Эйлера предлагается вычислять через другие промежуточные кинематические параметры — проекции x вектора конечного поворота x = —k = tg(p/4)e на связанные координатные оси по формулам
Л* = (1 — x2 + 2x)/(1 + x2), x2 = x • x = x2 + x2 + x2 (15)
или по эквивалентным им формулам
A* = (1 — x2 )/(1+ x2), Л* = 2x/(1+ x2), x2 = x? + x2 + x3 • (16)
В свою очередь, переменные x^ находятся по одному из выбранных алгоритмов численного интегрирования (с нулевыми начальными условиями) нового кватернионного кинематического уравнения типа Риккати (10) или (11), принимающего с учетом введенной новой переменной x вид (17) или (18):
4x^ = шх + x о шх — шх о x — x о шх о x, (17)
4x^ = шх — 2шх х x — x о шх о x. (18)
Соотношения (15) (или (16)), связывающие приращение А* кватерниона поворота А на шаге интегрирования с вектором поворота х, и дифференциальное уравнение (18) для вектора х гораздо проще соотношений (13), связывающих приращение А* с вектором поворота и дифференциального уравнения (14) для вектора ^ (уравнения Борца). Поэтому нами для построения алгоритмов ориентации БИНС предлагается использовать уравнение (18) и соотношения (15) (или (16)).
Приведем исследуемые новые алгоритмы вычисления приращений А* параметров Эйлера Аj на шаге интегрирования, построенные (с помощью метода последовательного приближения Пикара) с использованием формул Кэли и кватернионного кинематического уравнения типа Риккати [2-5].
2.1. Одношаговые алгоритмы третьего порядка точности
Алгоритм 1.
А* = (1 — х2 )/(1+ х2), А1 = 2х/(1 + х2), х2 = х? + х2 + х2,
х = а7 + 67* х 7 — 07 о 7* о 7, (19)
^ —1
7 = J 7* = J
Алгоритм 2. В этом алгоритме кватернионная переменная х с нулевой скалярной частью вычисляется по формуле
х = «7 + 67* х 7 + ¿(7*27* + 727), (20)
72 = 7 • 7 = 72 + 722 + 7э, 7*2 = 7* • 7* = 7*2 + 72*2 + 7з*2-Остальные соотношения алгоритма 1 не изменяются.
2.2. Двухшаговый алгоритм третьего порядка точности
А* = (1 — х2)/(1+ х2), А* = 2х/(1 + х2), х2 = х1 + х2 + х2,
11111111111
х = а(7 + 7 ) + в7 х 7 — £7 о 7 о 7 , (21)
7' = J 7'' = J ш(£)
—1 —1 + й/2
2.3. Двухшаговый алгоритм четвертого порядка точности
Первый вариант:
А* = (1 — х2)/(1+ х2), А* = 2х/(1 + х2), х2 = х1 + х2 + х2, х = а(У + 7 ') + еУ х 7" + 6(7 2 У + 7 '2 V) = (а + 67'2)У + (а + 67"2+ еУ х У'. (22) Второй вариант:
А* = (1 — х2 )/(1+ х2), А* = 2х/(1 + х2), х2 = х1 + х2 + х2, / » / // /2 / //2 // / // //2 / / // '2 // х = а(7 + 7 ) + е7 х 7 + 6(7 7 + 7 7 ) + 6[(7 ■ 7 + 7 )7 — (7 ■ 7 + 7 )7 ] =
' 2 ' '' ''2м 'г / ''2 ' '' '2м 11111
= [а + 6(7 2 + 7 ■ 7 + 7 2)]7 + [а + 6(7 2 — 7 ■ 7 — 7 2)]7 + е7 х 7 . (23)
В этих алгоритмах а, 6, с, е, I — числовые коэффициенты.
Замечание 1. Термины «одношаговый», «двухшаговый» и «четырехшаговый» алгоритм означают, что шаг интегрирования равен соответственно дискретности, удвоенной дискретности и учетверенной дискретности выдачи интегральной первичной информации гироскопов (измерителей абсолютной угловой скорости движущегося объекта).
3. АНАЛИЗ РЕЗУЛЬТАТОВ МОДЕЛИРОВАНИЯ
С применением математического моделирования работы БИНС исследованы известные (приведенные в [1]) алгоритмы определения инерциальной ориентации движущегося объекта, использующие интегральную первичную информацию о вращательном движении объекта: метод средней скорости [1, формулы (7.18)], имеющий второй порядок точности; одношаговый [1, формулы (7.20)] и двух-шаговый [1, формулы (7.22)] алгоритмы третьего порядка точности; двухшаговый алгоритм четвертого порядка точности [1, с. 355] и четырехшаговый алгоритм четвертого порядка точности А. П. Панова [9, с. 157, 158]:
А* = сов(р/2), Ак = р-1рк в1п(р/2) (к = 1, 2, 3), р = (р2 + р2 + р2)1/2, (24)
4
р = ^ 7« + 45[Г(1) + Г(2) ][7 (3) + 7(4) ] + 32[Г(1) 7(2) + Г(3) 7(4)], (25)
1=1
р = (р1 ,р2,р3), 7(г) = (7(г),7?г),73г)) (г = 1, 2, 3,4) — векторы-столбцы,
/ 0 (i) -y3 (i) y2 \
r(i) = (i) y3 0 (i -y1 ) ,
V (i) -y2 (i) y1 0 /
tn-i +h/4 tn -1 +h/2
тк1} = J U (t) dt, 7(2) = J uk
tn-1 tn-1+h/4
tn-l+3h/4 tn-1 +h
Yk3) = J uk(t) dt, 7k4) = J U(t)dt (k = 1,2,3), (26)
tn-i+h/2 tn-i +3h/4
получающийся (используются наши обозначения) из его алгоритма [9, формулы (3.3.38), (3.3.39), (3.3.45)] шестого порядка точности вычисления вектора ^ (опечатка: вместо (32/55) в формулах (3.3.45) из [9] должно быть (32/45)), а также новые алгоритмы ориентации, приведенные в данной статье: одношаговые алгоритмы третьего порядка точности (19), (20) и двухшаговый алгоритм четвертого порядка точности (22).
Отметим, что известный одношаговый алгоритм третьего порядка точности [1, формулы (7.20)] был приведен ранее в работах [13,14] (он также приводится в [12]), а двухшаговый алгоритм третьего порядка точности [1, формулы (7.22)] (см. также [12]) и двухшаговый алгоритм четвертого порядка точности были получены в [1, с. 355] на основе двухшагового алгоритма вычисления координат вектора ориентации ^ = ре четвертого порядка точности, приведенного в [9], при учете связей координат этого вектора с параметрами Эйлера. Эти алгоритмы ориентации отличаются простотой и хорошей точностью. Метод средней скорости, приведенный и исследованный ранее в [15-17], относится к алгоритмам второго порядка точности, которые широко использовались в алгоритмах БИНС космических кораблей «Союз Т/ТМ/ТМА», «Прогресс М/М1» [12]. Аналогичные алгоритмы второго порядка точности использовались и для орбитальных станций «Мир», Международной космической станции (МКС) [12]. Четырехшаговый алгоритм четвертого порядка точности А. П. Панова реализован в алгоритмах БИНС, построенных на лазерных гироскопах.
Одной из отличительных черт проведенного исследования заключается в изучении точности (вычислительного дрейфа) алгоритмов ориентации при наличии высокочастотных (порядка 100 или 400 гц) колебаний основания с учетом частоты съема первичной интегральной информации 1000 или 2000 гц.
Исследование проведено с помощью объемного математического моделирования. При моделировании движения объекта предполагается, что он совершает по каждой из угловых переменных 7 (типа самолетных углов) отдельные гармонические колебания или совершает по каждой из угловых переменных композицию нескольких гармонических колебаний с разными амплитудами и частотами. Задание таких законов движения объекта позволяет наиболее полно выявить точностные возможности исследуемых алгоритмов (традиционно при исследовании алгоритмов ориентации движение объекта
задается в виде конического движения, т.е. моделируется колебательное движение объекта по двум, а не трем угловым переменным). Вычисления проводились с 32- или 64-разрядной сеткой (последняя использовалась при моделировании алгоритмов четвертого порядка точности). Шаг интегрирования к выбирался из интервала [0.0005 с ^ 0.1 с]. Время движения объекта (интегрирования) — 600 с (10 мин).
Моделирование поводилось как для отдельных низкочастотных угловых гармонических колебаний объекта (частоты колебаний /^ = /7 = 1 гц, / = 0.5 гц) с малыми = 1 град, $+ = 2 град, 7+ = 3 град) и большими = 15 град, $+ = 5 град, 7+ = 15 град) амплитудами (а также для других вариантов отдельных низкочастотных угловых гармонических колебаний объекта с большими амплитудами и частотами), так и для композиций высокочастотных и низкочастотных гармонических колебаний объекта, в частности для следующих параметров гармоник:
первая высокочастотная гармоника: частоты — /^ = 380 гц, / = 400 гц, /7 = 420 гц; амплитуды (одинаковые) - = $+ = 7+ = 2.5 угл. с или = $+ = 7+ = 25 угл. с;
вторая высокочастотная гармоника: частоты — /^ =60 гц, / = 80 гц, /7 = 100 гц; амплитуды —
= 1 угл. мин, $+ = 1.5 угл. мин, 7+ = 2 угл. мин или = 10 угл. мин, $+ = 15 угл. мин, 7+ = 20. угл. мин;
третья низкочастотная гармоника: частоты — /^ = / = /7 = 1 гц (одинаковые) или /^ = 1 гц, / = 2 гц, /7 = 3 гц; амплитуды — = 6 град, $+ = 8 град, 7+ = 10 град.
Отметим, что приведенные ниже погрешности алгоритмов численного интегрирования, использующих интегральную информацию о движении объекта, включают в себя погрешности вычисления приращений интегралов от проекций вектора абсолютной угловой скорости объекта (т. е. включают погрешности формирования входной интегральной информации алгоритмов ориентации), величины которых зависят от числа М точек разбиения интервалов вычисления приращений интегралов (М = 11, 21,41,...) и от задаваемых угловых движений объекта (эта погрешность существенно влияет на точность определения ориентации с помощью алгоритмов третьего и выше порядков точности).
Приведем основные выводы, касающиеся точностных характеристик алгоритмов ориентации для различных угловых движений объекта (для отдельных трехчастотных гармонических колебаний объекта по углам Эйлера - Крылова и для их композиций).
3.1. Отдельные трехчастотные гармонические колебания объекта по углам Эйлера - Крылова
3.1.1. Случай, когда частоты колебаний /^ = /7 = 1 гц, / = 0.5 гц, а амплитуды колебаний либо малые: = 1 град, $+ = 2 град, 7+ =3 град, либо большие: = 15 град, $+ = 5 град, 7+ = 15 град.
1. Анализ изменения во времени методических погрешностей определения углов ориентации объекта показывает, что погрешность вычисления каждого из углов ориентации имеет две составляющие: медленно изменяющуюся составляющую, монотонно возрастающую во времени, и быстро изменяющуюся периодическую составляющую, период изменения которой определяется периодом гармонических колебаний объекта (для отдельных низкочастотных угловых гармонических колебаний объекта этот период равен 2 с или 1 с).
2. Погрешность вычисления углов ориентации при совершении объектом угловых колебаний с большими амплитудами увеличивается в среднем на два порядка по сравнению с погрешностями для угловых колебаний с малыми амплитудами.
3. При увеличении шага интегрирования на порядок погрешности вычисления углов ориентации увеличиваются на 2 ^ 3 и более порядков; при совершении объектом угловых колебаний с большими амплитудами худшие результаты для шага интегрирования к = 0.1 с дает одношаговый алгоритм третьего порядка точности [1, формулы (7.20)], использующий информацию из предыдущего шага интегрирования, а наименьшую — новый двухшаговый алгоритм четвертого порядка точности (22) и четырехшаговый алгоритм четвертого порядка точности А. П. Панова (24), (25); причем эти алгоритмы для этого большого шага интегрирования в отличие от других исследуемых алгоритмов остаются работоспособными, имея по угловым переменным ^,$,7 накопленные методические погрешности, равные 2.92 ■ 10-2 град, 4.52 ■ 10-3 град, 1.63 ■ 10-2 град и 2.09 ■ 10-2 град, 3.93 ■ 10-3 град, 1.48 ■ 10-2 град соответственно.
Замечание 2. Отметим, что при сравнении потенциальной точности двухшагового и четырех-шагового алгоритмов четвертого порядка точности нужно сравнивать результаты интегрирования, полученные с помощью четырехшагового алгоритма для шага интегрирования, вдвое превосходящего шага интегрирования дифференциальных уравнений ориентации объекта с помощью двухшагового алгоритма, так как двухшаговый алгоритм может быть реализован с частотой, вдвое превосходящей частоту реализации четырехшагового алгоритма. Для шага интегрирования h = 0.2 с четырехшаговый алгоритм (24), (25) имеет по угловым переменным $,y накопленные методические погрешности, равные 3.20 ■ 10-1 град, 1.20 ■ 10-1 град, 4.72 ■ 10-1 град соответственно. Эти погрешности более чем на порядок превосходят погрешности нового двухшагового алгоритма (22) для шага интегрирования h = 0.1 с.
4. Для моделируемых низкочастотных угловых движений объекта с малыми амплитудами метод средней скорости [1, (7.18)] имеет хорошую точность, поэтому он, как известно [12], широко используется для определения ориентации космических аппаратов. При шаге интегрирования h = 0.001 с его погрешности составляют величины 9.73 ■ 10-7 град, 1.25 ■ 10-7 град и 5.74 ■ 10-7 град по переменным
$, y соответственно, а при шаге интегрирования h = 0.01 с эти погрешности составляют величины 9.75 ■ 10-5 град, 1.20 ■ 10-5 град и 5.97 ■ 10-4 град соответственно (при одинаковом числе M = 11 точек разбиения шага интегрирования, используемых при вычислении интегралов от проекций угловой скорости объекта), т. е. увеличиваются на два порядка. Новый одношаговый алгоритм третьего порядка точности (19) имеет в этих случаях погрешности, равные 3.30 ■ 10-8 град, 1.95 ■ 10-8 град, 9.99■ 10-8 град и 4.65 ■ 10-6 град, 2.58■ 10-6 град, 1.62 ■ 10-5 град по переменным $, y соответственно (при том же числе M = 11 точек разбиения шага интегрирования).
Для моделируемых низкочастотных угловых движений объекта с большими амплитудами погрешности метода средней скорости увеличиваются в среднем на 2 порядка и составляют (при шагах интегрирования, равных 0.001 с и 0.01 с) величины 1.17 ■ 10-4 град, 4.33 ■ 10-5 град, 1.68 ■ 10-4 град и 1.20 ■ 10-2 град, 4.50 ■ 10-3 град, 1.72 ■ 10-2 град по переменным $, y соответственно, в то время как новый одношаговый алгоритм третьего порядка точности (19) имеет при шаге интегрирования h = 0.01 с значительно меньшие погрешности, равные 2.49 ■ 10-4 град, 2.67 ■ 10-3 град, 9.12 ■ 10-4 град по переменным $, y соответственно.
5. Для моделируемых низкочастотных угловых движений объекта с большими амплитудами из всех одношаговых алгоритмов ориентации наименьшую методическую погрешность имеет новый од-ношаговый алгоритм третьего порядка точности (19). При шаге интегрирования h = 0.001 с его погрешности составляют величины 3.53 ■ 10-7 град, 2.73 ■ 10-6 град и 4.12 ■ 10-7 град по переменным $, y соответственно, а при шаге интегрирования h = 0.002 с эти погрешности составляют величины 9.89 ■ 10-7 град, 2.14 ■ 10-5 град и 6.03 ■ 10-6 град соответственно, т.е. при увеличении шага интегрирования в 2 раза погрешность по переменной ^ увеличивается примерно в 3 раза, по переменной $ — в 8 раз, а по переменной y — в 15 раз. При шаге интегрирования h = 0.001 с погрешности известного одношагового алгоритма третьего порядка точности [1, (7.20)] составляют величины 6.29 ■ 10-7 град, 8.03 ■ 10-6 град и 1.26 ■ 10-6 град по переменным $, y соответственно, а при шаге интегрирования h = 0.002 с эти погрешности составляют величины 3.57 ■ 10-6 град, 6.42 ■ 10-5 град и 1.30 ■ 10-5 град соответственно, т.е. при увеличении шага интегрирования в 2 раза погрешность по переменной ^ увеличивается примерно в 6 раз, по переменной $ — в 8 раз, а по переменной y — в 10 раз. Число точек разбиения шага интегрирования, используемых при вычислении интегралов от проекций угловой скорости объекта, для шага интегрирования h = 0.001 бралось равным 21, а для шага интегрирования h = 0.002 бралось равным 41.
6. Из всех многошаговых алгоритмов ориентации наименьшую методическую погрешность имеет новый двухшаговый алгоритм четвертого порядка точности (22) и четырехшаговый алгоритм четвертого порядка точности А. П. Панова (24), (25). При шаге интегрирования h = 0.001 с погрешности двухшагового алгоритма (22) составляют величины 4.13 ■ 10-7 град, 1.47 ■ 10-7 град и 5.40 ■ 10-7 град по переменным $, y соответственно, а при шаге интегрирования h = 0.002 с эти погрешности составляют величины 1.66 ■ 10-6 град, 5.87 ■ 10-7 град и 2.16 ■ 10-6 град соответственно (при одинаковом числе M = 21 точек разбиения шага интегрирования, используемых при вычислении интегралов от проекций угловой скорости объекта). При шаге интегрирования h = 0.002 с погрешности четырехша-
гового алгоритма (24), (25) составляют величины 4.16 ■ 10-7 град, 1.47■ 10-7 град и 5.42 ■ 10-7 град по переменным 7 соответственно, а при шаге интегрирования к = 0.004 с эти погрешности составляют величины 1.70■ 10-6 град, 5.95-10-7 град и 2.19-10-6 град соответственно (при одинаковом числе М = 41 точек разбиения шага интегрирования, используемых при вычислении интегралов от проекций угловой скорости объекта). Погрешности нового двухшагового алгоритма (22) и четырехшагового алгоритма четвертого порядка точности А. П. Панова (24), (25) для этого случая движения объекта оказались примерно равными (с незначительным преимуществом нового алгоритма, см. замечание 2).
Погрешности известного алгоритма четвертого порядка точности [1, с. 355] при шаге интегрирования к = 0.001 с составляют величины 6.85 ■ 10-7 град, 2.64 ■ 10-7 град и 8.07 ■ 10-7 град по переменным 7 соответственно, а при шаге интегрирования к = 0.002 с эти погрешности составляют величины 4.83 ■ 10-6 град, 1.63 ■ 10-6 град и 5.31 ■ 10-6 град соответственно (при одинаковом числе М = 21 точек разбиения шага интегрирования, используемых при вычислении интегралов от проекций угловой скорости объекта).
Замечание 3. При реализации четырехшагового алгоритма четвертого порядка точности А. П. Панова (24), (25) может происходить потеря точности при переходе от вектора ориентации, вычисляемого по этому алгоритму, к параметрам Родрига - Гамильтона (Эйлера). Точные формулы этого перехода содержат такие трансцендентные функции, как синус, косинус, квадратный корень. Аппроксимация этих формул отрезками рядов с удержанием лишь величин третьего порядка малости приводит в анализируемом случае к существенной потере точности определения ориентации. Погрешности алгоритма в этом случае при к = 0.002 с равны 1.22 ■ 10-4 град, 6.03 ■ 10-5 град и 2.24 ■ 10-4 град по переменным 7 соответственно, что больше погрешностей одно- и двухшаговых алгоритмов на
3 порядка. Этот переход с использованием стандартных подпрограмм вычисления трансцендентных функций существенно повышает точность определения ориентации для моделируемых слабочастотных и высокочастотных угловых движений объекта с большими амплитудами и делает ее, как уже отмечалось, сравнимой с точностью определения ориентации с помощью нового двухшагового алгоритма четвертого порядка точности (22). Однако алгоритм (22) имеет преимущество в смысле объема необходимых вычислений в 1.5 ^ 2 раза из-за его большей простоты. Кроме того, двухшаго-вый алгоритм (22) позволяет выполнять определение ориентации с большей частотой (по сравнению с четырехшаговым алгоритмом (24), (25)) и, следовательно, большей точностью.
Замечание 4. Число точек разбиения шага интегрирования, используемых при вычислении интегралов от проекций угловой скорости объекта, существенно влияет на точность определения ориентации объекта. Так, при увеличении этого числа в 2 раза погрешности нового алгоритма (22) при шаге интегрирования к = 0.002 си М = 41 составляют величины 4.22 ■ 10-7 град, 1.50 ■ 10-7 град и 5.46 ■ 10-7 град соответственно, а погрешности алгоритма (24), (25) при шаге интегрирования к = 0.004 си М = 81 составляют величины 4.85 ■ 10-7 град, 1.6 ■ 10-7 град и 5.89 ■ 10-7 град соответственно. Видно (см. п. 6), что точность как первого, так и второго алгоритма существенно повысилась. Отметим, что в этом случае при увеличении шага интегрирования в 2 раза погрешности как первого, так и второго алгоритма определения ориентации увеличились незначительно, что говорит о меньшей чувствительности этих алгоритма к шагу интегрирования (об их лучшей вычислительной устойчивости).
3.1.2. Отдельные трехчастотные гармонические колебания объекта по углам Эйлера - Крылова с различными (большими) амплитудами и частотами.
В табл. 1 приведены результаты моделирования алгоритмов четвертого порядка точности (нового (22) и известного [1, с. 355] двухшаговых алгоритмов, а также четырехшагового алгоритма А. П. Панова (24), (25)) для следующих вариантов отдельных трехчастотных гармонических колебаний объекта по углам Эйлера - Крылова с большими частотами Шф, ш7, и амплитудами $+, 7+ на интервале времени I = [0,600с]:
1) Шф = ш7 = 2п рад/с, = п рад/с; = 15 град, $+ = 5 град, 7+ = 15 град;
2) Шф = ш7 = 2п рад/с, = п рад/с; = 30 град, $+ = 10 град, 7+ = 30 град;
3) Шф = ш7 = 4п рад/с, = 2п рад/с; = 15 град, $+ = 5 град, 7+ = 15 град;
4) Шф = ш7 = 4п рад/с, = 2п рад/с; = 30 град, $+ = 10 град, 7+ = 30 град;
5) Шф = ш7 = 2п рад/с, = п рад/с; = 60 град, = 20 град, 7+ = 60 град;
6) Шф = ш7 = 2п рад/с, = п рад/с; = 120 град, = 40 град, 7+ = 120 град.
Таблица 1
Погрешности алгоритмов для гармонических колебаний объекта с различными амплитудами и частотами
Метод интегрирования К, с М Вар. Погрешности вычисления углов, град
Аф А^ А7
Новый двухшаговый алгоритм четвертого порядка (22) 0.001 21 1 4.135 ■ Ю-7 1.468 ■ 10-7 5.401 ■ 10-7
2 3.787 ■ 10-6 1.387 ■ 10-6 2.617 ■ 10-6
3 6.446 ■ 10-7 1.954 ■ 10-7 6.673 ■ 10-7
4 1.513 ■ 10-6 7.630 ■ 10-7 1.551 ■ 10-6
5 6.886 ■ 10-5 3.835 ■ 10-6 1.231 ■ 10-5
6 1.826 ■ 10-3 5.421 ■ 10-4 7.067 ■ 10-4
0.002 21 1 1.657 ■ 10-6 5.873 ■ 10-7 2.161 ■ 10-6
0.01 21 1 4.390 ■ 10-5 1.494 ■ 10-5 5.501 ■ 10-5
2 4.156 ■ 10-4 1.342 ■ 10-4 2.521 ■ 10-4
3 6.409 ■ 10-5 1.898 ■ 10-5 6.541 ■ 10-5
4 1.498 ■ 10-4 8.944 ■ 10-5 1.432 ■ 10-4
0.1 21 1 2.922 ■ 10-2 4.522 ■ 10-3 1.635 ■ 10-2
Четырехшаговый алгоритм А. П. Панова четвертого порядка (24), (25) 0.002 41 1 4.164 ■ 10-7 1.473 ■ 10-7 5.420 ■ 10-7
2 3.814 ■ 10-6 1.387 ■ 10-6 2.617 ■ 10-6
3 6.443 ■ 10-7 1.947 ■ 10-7 6.669 ■ 10-7
4 1.511 ■ 10-6 7.717 ■ 10-7 1.545 ■ 10-6
5 6.944 ■ 10-5 4.544 ■ 10-6 1.289 ■ 10-5
6 1.885 ■ 10-3 6.009 ■ 10-4 7.832 ■ 10-4
0.004 41 1 1.704 ■ 10-6 5.953 ■ 10-7 2.193 ■ 10-6
0.02 41 1 7.352 ■ 10-5 1.996 ■ 10-5 7.475 ■ 10-5
2 6.922 ■ 10-4 1.378 ■ 10-4 2.521 ■ 10-4
3 6.358 ■ 10-5 2.355 ■ 10-5 6.482 ■ 10-5
4 1.657 ■ 10-4 1.756 ■ 10-4 1.145 ■ 10-4
0.2 41 1 3.190 ■ 10-1 1.193 ■ 10-1 4.721 ■ 10-1
Двухшаговый алгоритм четвертого порядка [1, с. 355] 0.001 21 1 6.850 ■ 10-7 2.640 ■ 10-7 8.072 ■ 10-7
2 1.958 ■ 10-5 5.941 ■ 10-6 1.829 ■ 10-5
3 1.072 ■ 10-5 3.796 ■ 10-6 1.080 ■ 10-5
4 3.558 ■ 10-4 1.217 ■ 10-4 3.560 ■ 10-4
5 5.528 ■ 10-4 1.854 ■ 10-4 4.956 ■ 10-4
6 1.376 ■ 10-2 6.380 ■ 10-3 1.240 ■ 10-2
0.002 21 1 4.827 ■ 10-6 1.631 ■ 10-6 5.311 ■ 10-6
0.01 21 1 5.397 ■ 10-4 1.752 ■ 10-4 5.501 ■ 10-4
2 1.673 ■ 10-2 5.644 ■ 10-3 1.656 ■ 10-2
3 1.121 ■ 10-2 3.774 ■ 10-3 1.122 ■ 10-2
4 3.613 ■ 10-1 1.234 ■ 10-1 3.614 ■ 10-1
Примечание. К — шаг интегрирования, с; М — число точек вычисления приращения интеграла.
Из анализа результатов моделирования, приведенных в табл. 1, следует, что для моделируемых угловых движений объекта наименьшие сопоставимые между собой погрешности имеют новый двух-шаговый алгоритм четвертого порядка точности (22) и известный четырехшаговый алгоритм четвертого порядка точности А. П. Панова (24), (25). Двухшаговый алгоритм (22) имеет небольшое преимущество по точности перед четырехшаговым алгоритмом (24), (25) (см., однако, замечания 2 и 3).
3.2. Композиции высокочастотных и низкочастотных угловых гармонических колебаний объекта по каждому из углов Эйлера - Крылова
3.2.1. Первая композиция высокочастотных и низкочастотных угловых гармонических колебаний объекта по каждой из угловых переменных, состоящая из первой высокочастотной гармоники (частоты — /^ = 380 гц, / = 400 гц, /7 = 420 гц; амплитуды — = §+ = 7+ = 2.5 угл. сек (одинаковые)), второй высокочастотной гармоники (частоты — /^ = 60 гц, / = 80 гц, /7 = 100 гц; амплитуды — = 1 угл. мин, §+ = 2 угл. мин, 7+ = 2 угл. мин) и третьей низкочастотной гармоники (частоты — /^ = / = / = 1 гц (одинаковые); амплитуды — = 6 град, §+ = 8 град, 7+ = 10 град); время движения Т = 600 с (табл. 2).
Таблица 2
Первая композиция высокочастотных и низкочастотных колебаний объекта
Метод интегрирования К, с М f, гц Погрешности вычисления углов, град
Аф А^ А7
Четырехшаговый алгоритм А. П. Панова четвертого порядка (24), (25) 0.002 41 2000 8.326 ■ 10-5 2.799 ■ 10-6 3.693 ■ 10-6
0.004 81 1000 1.654 ■ 10-3 1.416 ■ 10-5 5.031 ■ 10-6
Метод средней скорости (второго порядка) (7.18) [1] 0.0005 11 2000 2.461 ■ 10-6 2.98 ■ 10-6 4.102 ■ 10-6
0.001 21 1000 7.672 ■ 10-6 4.798 ■ 10-6 5.386 ■ 10-6
Одношаговый алгоритм третьего порядка (7.20) [1] 0.0005 11 2000 1.478 ■ 10-6 2.438 ■ 10-6 3.807 ■ 10-6
0.001 21 1000 2.784 ■ 10-6 3.617 ■ 10-6 4.109 ■ 10-6
Двухшаговый алгоритм третьего порядка (7.22) [1] 0.001 21 2000 1.273 ■ 10-6 2.236 ■ 10-6 3.697 ■ 10-6
0.002 41 1000 1.232 ■ 10-3 1.139 ■ 10-5 4.209 ■ 10-6
Одношаговый алгоритм третьего порядка (19) 0.0005 11 2000 1.442 ■ 10-6 2.332 ■ 10-6 3.808 ■ 10-6
0.001 21 1000 2.006 ■ 10-6 2.75 ■ 10-6 4.003 ■ 10-6
Двухшаговый алгоритм четвертого порядка (22) 0.001 21 2000 1.315 ■ 10-6 2.272 ■ 10-6 3.755 ■ 10-6
0.002 41 1000 1.232 ■ 10-3 1.113 ■ 10-5 4.793 ■ 10-6
Примечание. К — шаг интегрирования, с; М — число точек вычисления приращения интеграла; f — частота съёма интегральной первичной информации об угловом движении объекта, гц.
Из анализа результатов моделирования, приведенных в табл. 2, можно сделать следующие выводы.
1. Из одношаговых алгоритмов наименьшие методические погрешности имеет (для шагов интегрирования Н = 0.0005 с и 0.001 с) алгоритм (19), построенный на основе кватернионного уравнения типа Риккати. Двухшаговые алгоритмы имеют практически одинаковые погрешности, причем при увеличении шага интегрирования в два раза (с 0.001 с до 0.002 с) погрешность по переменной ^ увеличивается на 3 порядка, по переменной § — на 1 порядок, а по переменной 7 увеличивается не значительно. Таким образом, для этих параметров угловых движений объекта имеет место вычислительная неустойчивость по переменной ^ (эта неустойчивость обусловлена неустойчивостью вычисления параметра Эйлера Л2).
2. Наименьшие (примерно равные) методические погрешности имеют при шаге интегрирования 0.001 с двухшаговые алгоритмы, которые составляют по переменным $ и 7 величины, равные 1.3 ■ 10-6 град, 2.3 ■ 10-6 град и 3.7■ 10-6 град соответственно (отметим, что методические погрешности одношаговых алгоритмов при шаге интегрирования 0.0005 с имеют несколько большие величины).
3. Четырехшаговый алгоритм имеет (по сравнению с одно- и двухшаговыми алгоритмами) большие методические погрешности (см. замечание 2). Для шага интегрирования 0.002 с, реализуемом при частоте съема информации 2000 гц, погрешность этого алгоритма по переменной ^ больше почти на 1.5 порядка.
4. Методические погрешности метода средней скорости сопоставимы с погрешностями одно- и двухшаговых алгоритмов (имеют одинаковый порядок, равный 10-6).
3.2.2. Вторая композиция высокочастотных и низкочастотных угловых гармонических колебаний объекта по каждой из угловых переменных (частоты высокочастотных гармоник те же, что и для первой композиции, однако их амплитуды больше на порядок; частоты низкочастотной гармоники больше частот для первой композиции: /^ = 1 гц, / = 2 гц, /7 = 3 гц, амплитуды те же самые); время движения Т = 600 с (табл. 3).
Таблица 3
Вторая композиция высокочастотных и низкочастотных колебаний объекта
Метод интегрирования К, с М f, гц Погрешности вычисления углов, град
Д^ Д7
Четырехшаговый алгоритм А. П. Панова четвертого порядка (24), (25) 0.002 41 2000 8.220 ■ 10-3 5.870 ■ 10-5 4.087 ■ 10-5
0.004 41 1000 1.647 ■ 10-1 5.087 ■ 10-4 2.151 ■ 10-4
Метод средней скорости (второго порядка) (7.18) [1] 0.0005 11 2000 1.235 ■ 10-4 3.605 ■ 10-4 6.899 ■ 10-5
0.001 21 1000 5.038 ■ 10-4 1.369 ■ 10-3 1.652 ■ 10-4
0.002 41 1000 3.692 ■ 10-1 5.326 ■ 10-3 9.006 ■ 10-4
Одношаговый алгоритм третьего порядка (7.20) [1] 0.0005 11 2000 1.867 ■ 10-5 2.666 ■ 10-5 3.803 ■ 10-5
0.001 21 2000 5.715 ■ 10-5 4.378 ■ 10-5 4.869 ■ 10-5
0.002 41 1000 5.099 ■ 10-1 2.345 ■ 10-4 1.983 ■ 10-4
Двухшаговый алгоритм третьего порядка (7.22) [1] 0.001 21 2000 2.04 ■ 10-5 3.429 ■ 10-5 4.976 ■ 10-5
0.002 41 1000 1.231 ■ 10-1 1.140 ■ 10-4 1.417 ■ 10-4
Одношаговый алгоритм третьего порядка (19) 0.0005 11 2000 1.866 ■ 10-5 2.606 ■ 10-5 3.828 ■ 10-5
0.001 21 1000 5.715 ■ 10-5 4.378 ■ 10-5 4.869 ■ 10-5
0.002 41 1000 5.099 ■ 10-1 2.798 ■ 10-4 1.223 ■ 10-4
Одношаговый алгоритм третьего порядка (20) 0.0005 11 2000 1.866 ■ 10-5 2.522 ■ 10-5 3.835 ■ 10-5
0.001 21 1000 5.692 ■ 10-5 3.964 ■ 10-5 4.946 ■ 10-5
Двухшаговый алгоритм четвертого порядка (22) 0.001 21 2000 1.560 ■ 10-5 2.340 ■ 10-5 3.890 ■ 10-5
0.002 21 1000 1.229 ■ 10-1 1.193 ■ 10-4 1.658 ■ 10-4
Примечание. К — шаг интегрирования, с; М — число точек вычисления приращения интеграла; f — частота съёма интегральной первичной информации об угловом движении объекта, гц.
Из анализа результатов моделирования, приведенных в табл. 3, можно сделать следующие выводы. 1. Методические погрешности одношаговых алгоритмов (при Н = 0.0005 си Н = 0.001 с), а также двухшаговых алгоритмов (при Н = 0.001 с) увеличились на порядок по сравнению с погрешностя-
ми для первой композиции. При к = 0.002 с методические погрешности двухшаговых алгоритмов увеличились по переменной 7 в 30 раз, по переменной $ - на порядок, а по переменной ^ — на 2 порядка.
2. Наименьшие методические погрешности имеет новый двухшаговый алгоритм (22). При шаге интегрирования 0.001 с его погрешности составляют по переменным $ и 7 величины, равные 1.56 ■ 10-5 град, 2.34 ■ 10-5 град и 3.89 ■ 10-5 град соответственно (отметим, что методические погрешности одношаговых алгоритмов при шаге интегрирования 0.0005 с имеют несколько большие величины).
3. Четырехшаговый алгоритм (24), (25) имеет (для шага интегрирования 0.002 с, реализуемом при частоте съема информации 2000 гц) по сравнению с новым двухшаговым алгоритмом (22) (для шага интегрирования 0.001 с, реализуемом при той же частоте съема информации 2000 гц) большие методические погрешности (по переменной большие на 2 порядка) (см. замечание 2).
4. Для рассматриваемых параметров угловых движений объекта при к = 0.004 с для четырех-шагового алгоритма и при к = 0.002 с для всех других алгоритмов имеет место вычислительная неустойчивость по переменной ^ (при увеличении шага интегрирования с 0.002 с до 0.004 с для четырехшагового алгоритма и с 0.001 с до 0.002 с для всех других алгоритмов погрешность по этой переменной увеличивается на 3-4 порядка).
ЗАКЛЮЧЕНИЕ
По результатам проведенного математического моделирования можно сделать следующие основные выводы.
1. Из всех рассмотренных алгоритмов наименьшие методические погрешности имеет новый двухшаговый алгоритм четвертого порядка точности (22), построенный на основе нового кватернионного дифференциального кинематического уравнения типа Риккати. Этот алгоритм обладает лучшей вычислительной устойчивостью и имеет простую структуру. При шаге интегрирования к = 0.01 с его погрешности для моделируемых низкочастотных угловых движений объекта с большими амплитудами составляют 1.29-10-5град, 3.93-10-6 град и 1.45-10-5 град по переменным $ и 7 соответственно, что на 1 ^ 2 порядков меньше (для этого же шага интегрирования), чем погрешности всех других рассмотренных одно- и двухшаговых алгоритмов (см. п. 3.1.1). При шаге интегрирования к = 0.001 с погрешности двухшагового алгоритма (22) для этих угловых движений объекта составляют величины 4.13 ■ 10-7 град, 1.47 ■ 10-7 град и 5.40 ■ 10-7 град по переменным $ и 7 соответственно, а при шаге интегрирования к = 0.002 с эти погрешности составляют величины 1.66 ■ 10-6 град, 5.87 ■ 10-7 град и 2.16 ■ 10-6 град соответственно.
2. Из всех рассмотренных одношаговых алгоритмов ориентации наименьшую методическую погрешность имеет новый одношаговый алгоритм третьего порядка точности (19), построенный на основе кватернионного уравнения типа Риккати. Известные одно- и двухшаговый алгоритмы третьего порядка точности [1, формулы (7.20), (7.22)], а также новый двухшаговый алгоритм третьего порядка точности (21), имеют хорошие (сопоставимые между собой) точностные характеристики. Эти алгоритмы для большинства рассмотренных угловых движений объекта не намного уступают по точности алгоритму четвертого порядка точности (22) и, следовательно, с успехом могут быть использованы для решения ряда задач определения ориентации движущихся объектов.
3. Известный четырехшаговый алгоритм четвертого порядка точности А. П. Панова (24), (25) имеет точность, сопоставимую с точностью нового двухшагового алгоритма четвертого порядка точности (22), если фигурирующие в этом алгоритме трансцендентные функции (синус, косинус, квадратный корень) вычисляются по высокоточным алгоритмам (реализуемым, например, в стандартных подпрограммах вычисления этих функций). Однако новый двухшаговый алгоритм четвертого порядка точности имеет преимущество в смысле объема необходимых вычислений в 1.5 ^ 2 раза в силу его большей простоты. Кроме того, новый двухшаговый алгоритм позволяет выполнять определение ориентации с большей частотой (по сравнению с четырехшаговым алгоритмом) и, следовательно, с большей точностью.
4. Вычислительные дрейфы всех алгоритмов для высокочастотных колебаний объекта с малыми амплитудами (наиболее интересный для практики случай) имеют порядок 10-6 град.
5. Наличие высокочастотных составляющих в угловых колебаниях объекта может приводить к потере вычислительной устойчивости алгоритмов ориентации при неверном выборе шага интегрирования. Так, для второй из рассмотренных композиций высокочастотных и низкочастотных угловых гармонических колебаний объекта потеря вычислительной устойчивости алгоритмов ориентации наблюдалась уже при шаге интегрирования, равном 0.002 с.
Отметим, что выбор того или иного алгоритма определения ориентации движущегося объекта в конечном счете определяется классом движущихся объектов (параметрами углового движения объекта, временем его движения), требуемой точностью решения задачи определения ориентации и характеристиками используемого вычислителя.
Библиографический список
1. Челноков Ю. Н. Кватернионные и бикватернион-ные модели и методы механики твердого тела и их приложения. Геометрия и кинематика движения. М. : Физматлит, 2006. 511 с.
2. Челноков Ю. Н., Переляев С. Е., Челнокова Л. А. Новые уравнения и алгоритмы ориентации БИНС в четырехмерных кососимметрических операторах // Системный анализ, управление и навигация : сб. тез. докл. 14-й междунар. науч. конф. М. : Изд-во МАИ-ПРИНТ, 2009. С. 35-36.
3. Челноков Ю. Н., Переляев С. Е., Челнокова Л. А. Дифференциальные кинематические уравнения вращательного движения твердого тела в четырехмерных кососимметрических операторах и новые алгоритмы ориентации БИНС // Проблемы критических ситуаций в точной механике и управлении : материалы Всерос. науч. конф. с междунар. участием. Саратов : ООО Издат. центр «Наука», 2013. С. 315-320.
4. Челноков Ю. Н., Переляев С. Е. Новые уравнения и алгоритмы ориентации и навигации БИНС в четырехмерных кососимметрических операторах // Сб. материалов XXI Санкт-Петербургской меж-дунар. конф. по интегрированным навигационным системам. СПб. : ОАО «Концерн "ЦНИИ "Электроприбор», 2014. С. 308-312.
5. Переляев С. Е., Челноков Ю. Н. Новые алгоритмы определения инерциальной ориентации объекта // Прикладная математика и механика. 2014. Т. 78, вып. 6. С. 778-789.
6. Гантмахер Ф. Р. Теория матриц. М. : Наука, 1967. 576 с.
7. Stuelpnagel J. On the parametrization of the three-dimensional rotation group // SIAM Review. 1964. Vol. 6, № 4. P. 422-429.
8. Переляев С. Е. О соответствии трехмерных и четырехмерных параметров группы трехмерных вращений // Изв. РАН. МТТ. 2009. № 2. С. 30-44.
9. Панов А. П. Математические основы теории инерциальной ориентации. Киев : Наук. думка, 1995. 279 с.
10. Bortz J. E. A new mathematical formulation for strapdown inertial navigation // IEEE Transaction on Aerospace and Electronic Systems. 1971. AES-7. № 1. P. 61-66.
11. Savage P. G. Strapdown inertial navigation integration algorithm design. Pt. 1 : Attitüde algorithms // J. Guidance, Control and Dynamics. Vol. 21, № 1. 1998. P. 19-28.
12. Бранец В. Н. Лекции по теории бесплатформенных инерциальных навигационных систем управления. М. : МФТИ, 2009. 304 с.
13. Эдвардс А. Бесплатформенные инерциальные навигационные системы // Вопросы ракетной техники. 1973. № 5. С. 50-57.
14. Бесараб П. Н. Определение параметров пространственной ориентации движущегося объекта // Журн. вычисл. матем. и матем. физ. 1974. Т. 14, № 1. С. 240-246.
15. Бранец В. Н., Шмыглевский И. П. Применение кватернионов в задачах ориентации твердого тела. М.: Наука, 1973. 320 с.
16. Челноков Ю. Н. Об определении ориентации объекта в параметрах Родрига - Гамильтона по его угловой скорости // Изв. АН СССР. МТТ. 1977. № 3. С. 11-20.
17. Плотников П. К., Челноков Ю. Н. Сравнительный анализ точности алгоритмов определения ориентации объекта в параметрах Родрига - Гамильтона и направляющих косинусах // Космические исследования. 1979. Т. 17, вып. 3. С. 371-377.
An Investigation of Algorithms for Estimating the Inertial Orientation of a Moving Object Yu. N. Chelnokov1, S. E. Perelyaev2, L. A. Chelnokova3
1 Chelnokov Yurii Nikolaevich, Saratov State University, 83, Astrakhanskaya st., Saratov, Russia, 410012; RAS Institute of Precision Mechanics and Control, 24, Rabochayast., Saratov, Russia, 410028, ChelnokovYuN@info.sgu.ru
2Perelyaev Sergei Egororovich, "AeroSpetsProekt"Ltd., 35A, Gagarina st., Zhukovskii, Moscow Oblast, Russia, 140180, Pers2030@yandex.ru
3Chelnokova Lyudmila Aleksandrovna, RAS Institute of Precision Mechanics and Control, 24, Rabochaya st., Saratov, Russia, 410028, ChelnokovaLA@mail.ru
The new and known strapdown INS algorithms for high-precision estimation of the orientation parameters of a moving object (Rodrigues - Hamilton (Euler) parameters) in the inertial frame are investigated. The new algorithms are based upon using the classical Hamilton rotation quaternion, quaternion with zero scalar part, which is correlated to the classical rotation quaternion via the quaternion equivalent of Cayley formula, and also the new quaternion differential equation for the inertial orientation of a moving object. The new algorithms are developed using the Picard successive approximation method. These algorithms use the integral raw information about absolute angular motion of an object as input data. It is demonstrated that the new algorithms are superior to the known algorithms of the same order regarding accuracy and complexity.
Key words: moving object, Rodrigues-Hamilton (Euler) parameters, inertial orientation, Hamilton quaternion, quaternion matrix, Cayley formula, four-dimensional skew-symmetric operator.
References
1. Chelnokov Yu. N. Quaternion and Biquaternion Models and Methods of Mechanics of Solid Bodies and its Applications. Geometry and Kinematics of Motion. Moscow, Fizmatlit, 2006. 511 p. (in Russian).
2. Chelnokov Yu. N., Perelyaev S. E., Chelnokova L. A. New SDINS Equations and Algorithms for Orientation with Four-Dimensional Skew-Symmetric Operators. Book of abstracts, 14-th Intern. Scientific Conf. "System Analysis, Control and Navigation". Moscow, MAI-PRINT, 2009, pp. 35-36 (in Russian).
3. Chelnokov Yu. N., Perelyaev S. E., Chelnokova L. A. Differential Kinematic Equations of the Angular Motion of a Solid in Four-dimensional Skew-symmetric Operators and the new Strap-down INS algorithms for orientation. Proc. of the Conf. "Problems of Critical Situations in Precision Mechanics and Control". Saratov, OJSC "Nauka" Publ. Center, 2013, pp. 315-320 (in Russian).
4. Chelnokov Yu. N., Perelyaev S. E. New Equations and Algorithms of Orientation and Navigations for Strapdown INS with Four-dimensional Skew-symmetric Operators. Proc. of the XXI Saint-Petersburg Intern. Conf. on Integrated Navigation Systems. Saint-Petersburg, State Research Center of the Russian Federation, Concern CSRI Elektro-pribor, 2014, pp. 308-312 (in Russian).
5. Perelyaev S. E., Chelnokov Yu. N. New Algorithms for Estimating the Orientation of an Object. J. Ap-pl. Math. Mech., 2014, vol. 78, iss. 6, pp. 778-789.
6. Gantmakher F. R. Theory of Matrices. Moscow, Nauka, 1967, 576 p. (in Russian).
7. Stuelpnagel J. On the parametrization of the three-dimensional rotation group. SIAM Review, 1964, vol. 6, no. 4, pp. 422-430.
8. Perelyaev S. E. On the correspondence between the three- and four-dimensional parameters of the three-dimensional rotation group. Mech. Solids, 2009, vol. 44, iss. 2, pp. 204-213.
9. Panov A. P. Mathematical Background of Inertial Orientation Theory. Kiev, Naukova Dumka, 1995, 279 p. (in Russian).
10. Bortz J. E. A new mathematical formulation for strapdown inertial navigation. IEEE Trans. On Aerospace and Electronic Systems, 1971, vol. AES-7, no. 1, pp. 61-66.
11. Savage P. O. Strapdown inertial navigation integration algorithm design. Pt. 1: Attitude algorithms. J. Guidance, Control and Dynamics, 1998, vol. 21, no. 1, pp. 19-28.
12. Branetz V. N. Lectures on the Theory of Inertial Navigation Control Systems. Moscow, MFTI, 2009, 304 p. (in Russian).
13. Edwards A. Strapdown Inertial Navigation Systems. Rocket Technology, 1973, no. 5, pp. 50-57.
14. Besarab P. N. Estimation of the Orientation Parameters of a Moving Object. U.S.S.R. Comput. Math. Math. Phys., 1974, vol. 14, no 1, p. 242-248. DOI: 10.1016/0041-5553(74)90156-6.
15. Branetz V. N, Shmyglevsky I. P. Application of Quaternions in Problems of Orientation of a Rigid Body. Moscow, Nauka, 1973, 320 p. (in Russian).
16. Chelnokov Yu. N. On Estimating the Orientation of an Object in Rodrigues - Hamilton Parameters Using its Angular Velocity. Izv. Akad. Nauk. Mekh. Tverd. Tela, 1977, no. 3, pp. 11-20 (in Russian).
17. Plotnikov P. K., Chelnokov Yu. N. Comparative Analysis of Accuracy of Algorithms for Estimating the Orientation of an Object in Rodrigues-Hamilton Parameters and Direction Cosines. Cosmic Research, 1979, vol. 17, iss. 3, pp. 371-377.
