Моделирование вибрационной обстановки на борту космического аппарата с оценкой кинематической погрешности определения его угловой ориентации Текст научной статьи по специальности «Физика»

НАВИГАЦИОННЫЕ И ГИРОСКОПИЧЕСКИЕ!

СИСТЕМЫ |

УДК 621.396.988.6:629.19

Е. С. Л о б у с о в, А. Н. Ч у л и н

МОДЕЛИРОВАНИЕ ВИБРАЦИОННОЙ ОБСТАНОВКИ НА БОРТУ КОСМИЧЕСКОГО АППАРАТА С ОЦЕНКОЙ КИНЕМАТИЧЕСКОЙ ПОГРЕШНОСТИ ОПРЕДЕЛЕНИЯ ЕГО УГЛОВОЙ ОРИЕНТАЦИИ

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

E-mail: alexey.chulin@gmail.com

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

Рассмотрим задачу оценки точности алгоритма определения ориентации космического аппарата (КА) бесплатформенной инерциальной навигационной системой (БИНС) с учетом погрешностей вычислений, возникающих под влиянием вибрационной обстановки. Как показано в работах Л.Е.Гудмана и А.Р.Робинсона [1], а также В.Н.Бранца и И.П. Шмыглевского [2], при численном интегрировании кинематических уравнений углового движения в случае сложного пространственного движения возможно возникновение низкочастотной вычислительной погрешности БИНС, которая в работе [2] названа кинематической. Численные оценки этой погрешности для случая высокочастотных гармонических вибраций при некоторых вариантах вычислительных алгоритмов БИНС приведены в работах П. Дж. Саважа [3], Дж. Дж. Марка и Д.А. Тазартеса [4], Ю.А. Литмановича и др. [5].

Для минимизации этой погрешности Дж. Э. Борцем [6] предложена гибридная (аналого-цифровая) схема обработки измерений БИНС. Подобная (чисто цифровая) схема, где аналоговая часть заменена дискретной с высокой частотой дискретности, предложена П. Дж. Саважем [3].

Аналогичная погрешность возникает и в гиростабилизаторах, как показано в работе А.Ю. Ишлинского [7], где указанная погрешность называется кинематическим уходом.

Кинематическая погрешность и скорость ее изменения зависят от характеристик вибрационной обстановки на конкретном КА, поэтому эти характеристики необходимо учитывать при ее оценке. В то же время оценка вибрационной обстановки требует измерения и обработки данных в широкой полосе частот, что предъявляет повышенные требования в том числе и к быстродействию вычислительных средств БИНС. Таким образом, проблемным вопросом является обоснованный выбор алгоритма обработки данных и его частоты дискретности; при этом должны удовлетворяться требования минимального использования вычислительных ресурсов. Несмотря на полученные общие выражения для оценки кинематичекой погрешности, при разработке математического обеспечения БИНС для вновь разрабатываемого КА приходится выполнять значительный объем исследовательских работ, которые целесообразно заменить единой технологией проектирования. Исходя из потребности создания такой технологии, сформулируем цели настоящей работы:

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

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

Исследование влияния вибраций на вычислительную погрешность БИНС проводится в интересах проектирования БИНС. Методом исследования служит математическое моделирование.

Влияние вибрационной обстановки на вычислительные ошибки БИНС. Одной из задач БИНС КА является определение ориентации его связанной системы координат (ССК) ОХУ2 относительно некоторой инерциальной системы координат (ИСК) ОХоУ0^о путем численного решения кинематических уравнений углового движения КА, которые можно представить в различных формах, в частности относительно вектора истинного поворота (вектора Эйлера)

" + 1 х * +

Ф sin ф

1 —

пч О X (в х ш) (1)

2(1 - cos $)J v ' w

(в = — вектор истинного поворота от ИСК к ССК; $ — его модуль; £ — единичный вектор оси вращения; Ш - вектор угловой скорости КА в проекциях на оси ССК) или относительно кватерниона ориентации (параметров Родрига-Гамильтона)

Л = ^А о ш (2)

2

^A — кватернион перехода от ИСК к ССК, связанный с вектором

$ - $

истинного поворота соотношением Л = cos —+ <- sin —

2 2

При этом определение ориентации в БИНС можно рассматривать как результат численного решения кинематических уравнений в форме (1) или (2) с использованием измеренного вектора ш.

Таким образом, задача синтеза алгоритма определения ориентации в основном состоит в выборе метода численного интегрирования уравнения (1) и частоты его дискретности в зависимости от требуемой точности знания ориентации, а также с учетом вычислительных возможностей аппаратуры БИНС.

Как показано в работе [2], одной из составляющих суммарной погрешности определения ориентации БИНС является кинематическая погрешность J, представляемая в форме

ф

j = л о [me х de*E] о Л,

(3)

где Ф — суммарный угол поворота; 89^ = 8шЕdr — вариация вектора

to

кажущегося поворота в*Е = j шdт (квазикоординат) в проекциях на

¿0

оси ССК; 6йЕ — вариация вектора ш в проекциях на оси ССК; d9*E — дифференциал вектора кажущегося поворота в проекциях на оси ССК; Л — кватернион, сопряженный Л.

Очевидно, что эта погрешность равна нулю при плоском вращении, т.е. когда 6в*Е и d9*E коллинеарны; если вращение носит пространственный характер, т.е. траектория вращения в квазикоординатах в*Е(£) представляет собой некоторую произвольную кривую С, то оценка "сверху" кинематической погрешности может быть определена по теореме Стокса [2]:

J II ^

Ö6E X de*E

C

Ö6E X de*E

7UC

= 2s,

(4)

где С' — дополнение кривой С до замкнутого контура, такое, что / 8в*Е х d#E = 0; 5 — площадь поверхности, охваченной контуром

C'

интегрирования C U C'.

о

t

t

Погрешность такого рода появляется, например, в случае конической прецессии твердого тела [2], который можно рассматривать как частный случай пространственных угловых колебаний: пусть вектор угловой скорости в проекциях на оси ССК представляется в виде

ш = a (i cos Ш + j sin Ш), (5)

где a, П — константы; i, j — орты ССК (проекция Ш на третий орт равна нулю).

В работе [2] показано, что в этом случае система смещается за период прецессии относительно исходного положения на

\\ц = (П у п, (6)

а средняя угловая скорость смещения (дрейфа)

a2

шдр = 2П, (7)

что совпадает с оценкой кинематической погрешности по выражению (5). Так, для прецессии с амплитудой a/П = 10-2 рад и круговой частотой П = 628 с-1 (что соответствует частоте П/2п = 100 Гц)

средняя угловая скорость смещения составит в наихудшем случае

a2

шдр = — = 0,0314рад/с, что недопустимо для большинства прак-

2П

тических приложений.

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

[1, 3, 6, 8].

Таким образом, при воздействии на вход БИНС угловой скорости как изменяющегося во времени вектора (что имеет место, в частности, при действии вибраций) может возникать кинематическая погрешность, имеющая характер дрейфа, для устранения которой некоммутативные составляющие в уравнениях кинематики (1) или (2) должны учитываться алгоритмами БИНС. Очевидно также, что для корректного воспроизведения некоммутативных составляющих при действии вибраций спектр угловой скорости должен лежать в полосе пропускания БИНС, а соотношение верхней граничной частоты спектра угловой скорости и частоты дискретности БИНС — удовлетворять условию теоремы Котельникова.

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

Для всех рассматриваемых вариантов примем, что интегрирование уравнений в БИНС осуществляется в конечных разностях [2, 9]. В этом случае на первом этапе обработки измеренный вектор угловой скорости в осях ССК ш преобразуется в оценку вектора кажущегося поворота

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

На втором этапе численно решаются уравнения кинематики одним из нижеперечисленных одношаговых или многошаговых методов.

В качестве одношаговых методов рассматриваются метод на основе разложения Пикара 2-го порядка (модифицированный метод Эйлера [2, 9]) и метод Рунге-Кутты 4-го порядка.

В качестве многошаговых методов рассматриваются неявный метод Адамса 2-го порядка (неявный метод трапеций); метод на основе разложения Пикара 3-го порядка [2, 9]; двухчастотный метод Са-важа [3].

Характеристика действующих вибраций. В качестве исходных данных рассматриваются два различных по частотным свойствам типа вибраций, характерных для пассивных и активных участков полета КА.

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

На активных участках полета КА основным источником вибрационных возмущений является двигательная установка (ДУ) и ее агрегаты. В этом случае возмущения, как правило, носят характер широкополосного случайного шума. На рис. 1 показан пример спектральных плотностей линейных вибрационных ускорений КА под действием ДУ по осям ССК Sax, Say, Saz как функций частоты f, оцененных по полетной телеметрии методом быстрого преобразования Фурье с осреднением по 100 выборкам. Как следует из рисунка, спектр вибраций существенно отличается от спектра белого шума и имеет резкие локальные максимумы.

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

(8)

Рис. 1. Пример спектральной плотности случайных вибрационных ускорений КА под влиянием работы ДУ (по одной из осей ССК)

Гармонические вибрации задаются в форме углов последовательных поворотов как одновременно действующие угловые колебания вокруг осей Y и Z с одинаковой частотой и разностью фаз, равной п/2 (что соответствует рассмотренному случаю конической прецессии):

ф^) = фт sin Ш;

9(t) = 9m cos Ш;

Y(t) = 0,

где ф, 9, y — углы последовательных поворотов (ф — угол рыскания, 9 — угол тангажа, y — угол крена), соответствующие последовательности поворотов вокруг осей Y—Z—X. Амплитуды фт, 9m и частота Q колебаний выбираются в соответствии с расчетными оценками для проектируемого КА.

Угловая скорость вычисляется из кинематических уравнений в углах последовательных поворотов ф, 9, Y:

Y + ф sin 9 ф cos 9 cos y + 9 sin y • (9)

9 cos y — ф cos 9 sin y

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

Шх

Ш = Шу =

Шг

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

я

ё(*п) = £ Хг (Ьп), (10)

Г=1

где Ьп — дискретное время; хг — вектор вибрационного ускорения г-го тона колебаний, представляемый в виде

Хг (Ьп) = 2 [Й1г (Ьп) ЙШ^Г 1п + Й2т (Ьп) СОБ РгЬп] , (11)

рг — несущая круговая частота г-го тона колебаний.

Модулирующие функции каждого из тонов й1г, й2г задаются в виде случайных непрерывных кусочно-линейных функций времени:

Ulr (tn) = Cr < Xir (Пг) + [X^ (Пг + 1) - Xiir (n

"tn Пг ^^^r

T

r

r

nr Tr tn

< (nr + 1)Tr;

(12)

f - - tn — щ. Tr

U2r (tn) = CA X2r (Пг ) + [X^2r (Пг + 1) — X^2r (Пг^ П

nrTr ^ tn < (nr + 1)Tr,

где nr = 0,1, 2,...; Cr — константа (весовой коэффициент); Xjr(nr) =

■ ix

i-ijr ( n J

= [ A1jr A2jr A3jr ]т (Ajjr(n) ^ N(0, 1), i = 1, 2, 3) — векторный дискретный белый гауссов шум; Tr — шаг отсчетов r-й модулирующей функции.

Такт дискретности модели T0 выбирается из условия теоремы Ко-тельникова

T0 < min Tr/2.

r

В этом случае спектральная плотность каждого из тонов имеет вид

S,,, (v) = (CrTr)2 p^^TTß (13)

[(v — Vr) Tr/2]

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

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

S££(V) = £ Sxrxr (v) = £ (CrTr)2 sin4 (v— /2. (14)

r=1 r=1 [(v — Vr) Tr/2\

Используя данный метод и варьируя число тонов Л и их параметры рг, Сг, Тг, можно получить семейство псевдослучайных функций вида (11), каждая из которых будет соответствовать спектру вида (13). При этом сами параметры Л, рг, Сг, Тг следует выбирать таким образом, чтобы спектр модели был максимально приближен к спектру реальных вибраций, а именно: Л равно числу резонансных пиков реального спектра, уг — резонансной круговой частоте г-го пика, параметр Тг должен задаваться как величина, обратная ширине г-го пика, а Сг — подбираться так, чтобы мощность модели в окрестности рг была бы близка к мощности реального спектра в той же окрестности.

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

Моделирование. Оценка кинематической погрешности БИНС выполнена путем математического моделирования работы БИНС при воздействии на объект угловых колебаний, соответствующих по своим спектральным свойствам реальным вибрационным возмущениям, действующим на КА в полете. Блок-схема моделирования показана на рис. 2. В результате моделирования оценивается вычислительная ошибка БИНС в виде углов малых поворотов Д0, ДО, Д7 вокруг осей X, У, ^ соответственно. Эти углы определяются по рассогласованию между кватернионом ориентации КА Л, вычисляемым эталонной моделью кинематики углового движения КА, и оценкой кватерниона

Рис. 2. Блок-схема моделирования

шхдр = ———, где Т — интервал моделирования, как функции круго-

ориентации КА Л, вычисляемой моделью алгоритма интегрирования кинематических уравнений в БИНС, при одном и том же векторе угловой скорости КА Ш, поступающей на вход. Эталонная модель реализует интегрирование кинематических уравнений КА неявным методом трапеций с высокой частотой дискретности (6,4 кГц), что существенно превышает верхнюю граничную частоту рассматриваемого спектра вибраций. Вследствие этого вычислительная погрешность эталонной модели пренебрежимо мала. В модели алгоритма интегрирования кинематических уравнений в БИНС реализуется каждый из методов интегрирования, представленных ранее, с возможностью варьирования частоты дискретности от прогона к прогону.

Модель вибрационной обстановки имитирует описанные гармонические или случайные вибрации объекта. Гармонические вибрации задаются по двум взаимно ортогональным осям (в качестве таковых условно выбраны оси У и Z ССК) как гармонические функции времени, при этом разность фаз между ними принимается равной п/2. По результатам моделирования оценивается ошибка БИНС Д7 (Т) по третьей (X) оси ССК и средняя скорость дрейфа по этой оси

Дт(Т) Т

вой частоты вибраций П. Дополнительно оценивается относительная

П

ошибка БИНС £ = 2— |шждр | как функция отношения частоты вибра-

а2

ций 2пП к частоте дискретности (относительной частоты вибраций).

Случайные возмущения задаются по трем осям ССК при этом оцениваются ошибки и скорости дрейфа БИНС по осям ССК.

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

На рис. 3 показана зависимость погрешности определения ориентации от метода интегрирования при постоянной частоте квантования и варьируемой частоте вибраций для двух случаев. В первом случае (рис. 3, а) частота квантования более чем в 2 раза превышает частоту вибраций, т.е. выполняется условие теоремы Котельникова. Во втором случае (рис.3,б) условие теоремы Котельникова не выполняется. В обоих случаях амплитуда вибраций по каждой оси принимается постоянной и равной 1%.

На рис. 4 показана зависимость относительной погрешности определения ориентации от относительной частоты вибраций в логарифмическом масштабе.

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

Анализ результатов моделирования. Как следует из рис. 3-5, при воздействии вибрационных возмущений вокруг двух осей (У и Z)

Рис. 3. Зависимость погрешности определения ориентации БИНС от частоты гармонических вибраций при различных методах интегрирования и фиксированной частоте квантования, а также при выполнении (а) и невыполнении (б) условия теоремы Котельникова

выходная информация БИНС содержит низкочастотную ошибку по третьей оси (т.е. по каналу Х).

Сравнение результатов, приведенных на рис. 3, показывает, что при выполнении условия теоремы Котельникова погрешность существенно меньше, чем при его невыполнении; при этом в последнем случае погрешность близка к аналитической оценке "сверху" по формуле (8). Исходя из этого, делаем вывод, что выбирая частоту дискретности таким образом, чтобы она более чем в 2 раза превышала частоту возмущений, можно добиться снижения скорости дрейфа.

Как видно из рис. 3, а, при относительно низких частотах возмущения (значительно меньших, чем частота Котельникова /о/2, где /0 — частота дискретности БИНС) скорость дрейфа в целом снижается с повышением порядка метода, это объясняется тем, что метод более высокого порядка при прочих равных условиях дает лучшую аппроксимацию правой части дифференциальных уравнений кинематики.

Рис. 4. Зависимость относительной погрешности определения ориентации БИНС от относительной частоты гармонических вибраций (обозначения — см. рис. 3, а

6.00Е-05

5.00Е-05

4.00Е-05

9- 3.00Е-05

2.00Е-05

1.00Е-05

0.00E+0G

[

м *

* i 1 1 \

ft-\

V Чх

--- - г" —ч

100

200

300

400

Частота дискретности, Гц

Рис. 5. Зависимость погрешности определения ориентации БИНС от частоты дискретности при случайных вибрациях (обозначения — см. рис. 3, б)

В то же время при использовании неявного метода трапеций скорость дрейфа нарастает с ростом частоты возмущений относительно медленно по сравнению с применением явного метода даже более высокого (3-го) порядка (см. рис.3,а и 4). Кроме того, при использовании неявного метода трапеций график скорости дрейфа пересекает ось

абсцисс при значении частоты возмущений, приблизительно равным 0,23/о (рис. 3, а, 5), т.е. в этой точке скорость дрейфа будет равна нулю. Это свойство может быть использовано для компенсации дрейфа БИНС, возникающего вследствие вибраций на определенной частоте.

При невыполнении условий теоремы Котельникова (см. рис. 3, б), методы интегрирования разных порядков практически не различаются по точности. Исключение представляет лишь метод Рунге-Кутты 4-го порядка, для которого характерен существенно более низкий по сравнению с другими методами уровень дрейфа. Это объясняется тем, что при использовании метода Рунге-Кутты частота дискретности измерений превышает частоту дискретности интегрирования (для метода 4-го порядка — в 2 раза), что расширяет диапазон выполнения условия теоремы Котельникова.

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

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

СПИСОК ЛИТЕРАТУРЫ

1. Goodman L. E., Robinson A. R. Effects of finite rotations on gyroscope sensing devices // J. of Applied Mechanics. - 1958. - Vol. 25. - P. 210-213.

2. Б p а н е ц В. Н., Ш м ы г л е в с к и й И. П. Введение в теорию бесплатформенных инерциальных навигационных систем. - М.: Наука, 1992. - 280 с.

3. Savage P. G. Strapdown analytics. - Maple Plane, NM: Strapdown Associates, Inc., 2000. - 1532 c.

4. M a r k J. G., T a z a r t e s D. A. Tuning of coning algorithms to gyro data frequency response characteristics // J. of Guidance, Control, and Dynamics. - 2001. - Vol.24, no. 4. - P. 641-647.

5. Литманович Ю. А., Лесючевский В. М., Гусинский В. З. Обработка информации с использованием приращений кратных интегралов от сигнала: От бесплатформенных ИНС к другим системам реального времени // Гироскопия и навигация. - 2000. - № 2 (29). - C. 46-56.

6.BortzJ. E. A new concept in strapdown inertial navigation. NASA Technical Report. - Washington, D.C.: NASA, 1970. - 224 c.

7. Ишлинский А. Ю. Ориентация, гироскопы и инерциальная навигация. -М.: Наука, 1976. - 672 с.

8. T i 11 e r t o n D. H., W e s t o n J. L. Strapdown inertial navigation technology. -The Institution of Electrical Engineering, 2004. - 558 c.

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

10. А н у ч и н О. Н., К о м а р о в а И. Э., П о р ф и р ь е в Л. Ф. Бортовые системы навигации и ориентации искусственных спутников Земли. - СПб.: ГНЦ РФ ЦНИИ "Электроприбор", 2004. - 326 с.

Статья поступила в редакцию 21.06.2010

Евгений Сергеевич Лобусов — канд. техн. наук. доцент кафедры "Системы автоматического управления" МГТУ им. Н.Э. Баумана.

Ye.S. Lobusov — Ph. D. (Eng.), assoc. professor of "Systems of Automatic Control"

department of the Bauman Moscow State Technical University.

Алексей Николаевич Чулин — начальник сектора ФГУП "НПО им. С.А. Лавочкина".

A.N. Chulin — head of sector of the Federal State Unitary Enterprise "NPO imeni S.A.

Lavochkina"