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

УДК 629.78

МЕТОДИКА ВОССТАНОВЛЕНИЯ ОРИЕНТАЦИИ КОСМИЧЕСКОГО АППАРАТА ПРИ КОМПЛЕКСИРОВАНИИ МАГНИТОМЕТРИЧЕСКИХ И РАДИОНАВИГАЦИОННЫХ ИЗМЕРЕНИЙ

© 2007 И. В. Белоконов, А. В. Крамлих

Самарский государственный аэрокосмический университет

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

Введение

Методы решения задачи определения ориентации космических аппаратов (КА) по магнитометрическим измерениям изложены в работах [1-4]. Основным недостатком этих методов является использование модели движения КА, что, в свою очередь, затрудняет использование этих методов для определения ориентации в темпе поступления информации. В работах [5-8] описаны методы, базирующиеся на согласовании векторов в двух системах координат (СК), при этом минимальное количество векторов равно двум [7]. В качестве системы векторов наиболее часто используются векторы напряженности магнитного поля Земли (МПЗ), направления на звезды и Солнце [6, 8]. Компоненты этих векторов, заданные в удобной для КА системе координат (например, в орбитальной СК), определяются, исходя из имеющейся априорной информации. В частности, вектор напряженности МПЗ отыскивается с использованием модели МПЗ, а векторы направлений на звезды и Солнце находятся по каталогам и данным об эфемеридах. Компоненты векторов в связанной с КА системе координат измеряются с помощью трехкомпонентного магнитометра и специальной аппаратуры, определяющей направления на звезды и Солнце. При определении ориентации КА на основе согласования векторов в двух СК необходимо знание орбиты движения КА, что требует наличия на нем навигационного приемника (НП) или радиоконтроля орбиты.

В работе предлагается методика определения ориентации и динамики движения КА с использованием минимального состава измерительной аппаратуры, в качестве ко -торой используется многоканальный НП, принимающий сигналы от спутниковых радионавигационных систем (СРНС) ГЛО-НАСС и GPS, и магнитометр.

Постановка задачи определения ориентации КА

При постановке и решении задачи определения ориентации КА использованы правые ортогональные СК с центром, расположенным в центре масс:

- связанная СК (ССК) OX1Y1Z1 (ось OX1 - продольная ось);

- орбитальная СК (ОСК) OX2Y2 Z2 (ось

OZ2 направлена по радиусу-вектору КА, ось

OY2 направлена по вектору кинетического момента орбитального движения КА, ось OX2 дополняет систему до правой).

Положение СК OX1Y1 Z1 относительно СК OX2Y2Z2 задается с помощью кватерниона v = (v0 ,v1 ,v2 ,v3), имеющего единичную

норму: v^ + v22 + v2 + v^ = 1. Матрицу перехода от OX2Y2Z2 к OX2Y2Z2 обозначим MXX2.

Элементы этой матрицы выражаются через компоненты кватерниона v с помощью формул:

2222 mii = vo + vi -V2-v3;

m

12

m

1з

m

21

: 2 • (viV2 + V0V3 )>'

:2 ^(viv3- V0V2);

: 2 • (viv2- vov3);

m22 = V02 - V12 + V22 - V32 ;

: 2 •(v0vi + V2V3 );'

2 •(viv3 + V0V2);

: 2 • (v2V3 - V0V2 );

(1)

m

23

m

31

m

32

2222 m33 = V0 - V1 - V2 + V3 •

Задача определения ориентации КА рассматривается как задача нахождения кватерниона V.

При разработке алгоритмов решения задачи определения ориентации широко применяется подход, основанный на согласовании измерений различных векторов в двух СК, взаимная ориентации которых подлежит определению [7-9]. При решении задачи определения ориентации в качестве первого

вектора и1 будет взят вектор напряженности МПЗ, а в качестве второго вектора и2 -вектор положения антенны НП. Определение

вектора и2 принципиально возможно по анализу пространственного расположения видимых и невидимых навигационных спутников (НС).

Для отыскания кватерниона используется метод, описанный в [8]. Суть метода заключается в минимизации критерия, представляющего собой взвешенную с весами а, сумму квадратов разностей между значениями двух векторов, заданных в двух СК [9]:

2

^ Мхх, )= £ а,- (и 1 - мХ1Х2 • и 2) (и - мХ1Х2 • и 2),

,=1

(2)

где мХ1Х2 - матрица, описывающая связь ОСК и ССК, параметризованная с помощью

кватернионов; и, и2 - векторы в ССК и ОСК, соответственно (, = 1,2).

После отыскания кватерниона V проекции абсолютной угловой скорости Ш СК ОХ^1'£1 на ее собственные оси находятся с

помощью численного дифференцирования найденного кватерниона и кинематических уравнений

W = 2(v0V1 - V1V0 + V3V2 - V2V3 )-W2 = 2(v0V2 - V3V0 + V1V3 - V3V1 )-W3 = 2(V0V3 - V3V0 + V2V1 - V1V2 )

(3)

Решение задачи определения ориентации КА

Решение задачи определения ориентации КА разбивается на два этапа.

На первом этапе отыскивается вектор положения антенны НП в ОСК, и с этой целью анализируется пространственное положение НС систем ГЛОНАСС и GPS. Все НС разделяются на видимые и невидимые, которые, в свою очередь, разделяются на невидимые из-за затенения Землей и затененные конструкцией КА.

Для определения вектора положения антенны НП в ОСК предполагается, что заданы следующие исходные данные:

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

2. Навигационные данные, формируемые НП (массив номеров всех навигационных спутников, массив номеров видимых НС, массив номеров невидимых НС, геоцентрические координаты всех НС в СРНС, представленные в виде матрицы размером КНСх3).

3. Параметры движения центра масс (ПДЦМ) КА, получаемые от НП.

По имеющимся исходным данным в ОСК вычисляются единичные векторы, кол-линеарные векторам дальностей до видимых (В) и невидимых (НВ) НС, и из них формируются соответствующие матрицы НВ размером ЫНВ х 3 и ННВ размером ЫННВ х 3 (при этом исключаются из рассмотрения те НС, видимость которых отсутствует из-за затенения Землей).

Введем обозначения:

J,

НВ = \grad В1 grad В

= ^gradНВ1 gradНВ2 K gradHВNВ Г -

Н

где gradi = {х2і,у2і,г2і} - единичный вектор дальности до і-го НС в проекциях на оси ОСК.

Исходя из того, что ширина диаграммы направленности антенны составляет 180°, для видимых и невидимых НС выполняем следующие соотношения:

cos (u12 , gradв ) > О, (i = 1-NB );

cos(u12,gradнві )< a (j = 1-NнB)

где U12 ={x2 , У2 ,z2 1

единичный вектор антенны, записанный в проекциях на оси ОСК.

Так как

U1

=1 и

gradt = 1, то, пред-

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

X2iX2 + У2іУ2 + Z2iZ2 > 0- (i = I NВ )

0 0 = 1,NH3 )

X2jX2 + У2jy2 + Z2jZ2 <

(4)

Используя соотношения (4), описывающие геометрические связи между видимыми и невидимыми НС и вектором антенны, можно записать функционал вида:

Ф(

B

X2 , У2 , Z2 ) = і (x2iX2 + У2іУ2 + Z2iZ2 - 1)2

i =1

NHB ,

+ i(x2 jX2 + У2 У + Z2 jZ2 +

j=1

1)2,

(5)

который в дальнейшем используется для поиска координат антенны в ОСК.

Первое слагаемое функционала описывает связь проекции вектора антенны и2 с проекциями единичных векторов видимых

На втором этапе непосредственно решается задача определения ориентации и динамики КА.

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

2 | 2 . 2 . 2 1

элементов кватерниона: v0 + V! + v2 + Vз = 1. В работе [8] показано, что минимизация критерия (2) при условии нормировки для элементов кватерниона сводится к нахождению минимального собственного числа четырехмерной матрицы:

21 B = і -“1 —

' S ZT

Z

t

(б)

где 5=I((и; )- и'(и,') - и; (к }:

I = -(щ хЩ'); I = -1р;)Щ.

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

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

/к)

(к-1 Vk) > 0

(k = 1,N).

i=0

НС gradВ на оси ОСК, второе слагаемое

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

Решается задача отыскания минимума

функционала (5) по координатам х2,у2,г2 с учетом условия нормировки для координат

2 ■ 2 і 2 і антенны: х2 + у2 + г2 = 1.

После уточнения знака кватерниона определяются проекции абсолютной угловой

скорости W СК OX2Y2Z2 на ее собственные оси по соотношениям (3).

Описание модельной задачи

Моделирование задачи определения ориентации КА проводилось при следующих положениях:

1) орбита КА круговая, высота (h) 300 км и 1000 км, наклонение 63°;

2) количество НС равно 48, что соответствует общему количеству НС в СРНС ГЛОНАСС (при ее полном развертывании) и GPS;

0

3) антенна НП расположена по оси

OX1;

4) вектор напряженности МПЗ считается точно измеренным;

5) положение КА на орбите задается случайным образом по равновероятному закону (от 0° до 3б0°);

6) массив углов ориентации КА формируется случайным образом по равновероятному закону, углы ориентации изменяются от 0° до 3б0°.

Моделирование задачи определения ориентации и динамики движения КА проводилось в три этапа.

Этап 1. Моделирование СРНС ГЛО-НАСС и GPS (для простоты моделирования предполагалось, что в каждый момент времени положения ГЛОНАСС/GPS спутников «заморожено»). Моделирование движения КА и магнитометрических измерений. Моделирование магнитометрических измерений проводилось следующим образом. В ОСК по модели МПЗ в виде модели прямого диполя [10] рассчитывался вектор напряженности

МПЗ ( U2l ), а затем с использованием известной матрицы перехода MX1X2 в ССК рассчитывался «измеренный» вектор напряженности МПЗ (Ul). Вектор антенны НП (Ul2) в ССК согласно допущениям задавался векто-

ром с координатами {1,0,0}. Исключались НС, невидимые из-за затенения Землей.

Этап 2. Непосредственное отыскание вектора антенны НП в ОСК, основанное на отыскании минимума функционала (5) по

координатам х2,у2,22 с учетом условия нормировки для координат антенны:

2 I 2 ■ 2 1

Х2 + У 2 + ^2 = 1.

Этап 3. Определение ориентации КА по комплексированию магнитометрических и радионавигационных измерений.

Исследование эффективности на модельной задаче

Для исследования эффективности решения задачи определения ориентации КА при комплексировании магнитометрических и радионавигационных измерений была сформирована выборка решений объемом 100000 реализаций.

Для высот 300 и 1000 км построена плотность распределения ошибки положения

антенны Р(йа) (рис. 1 и 2). В качестве погрешности определения вектора положения антенны взят пространственный угол (да) между истинным и найденным вектором положения антенны НП.

Математическое ожидание ошибки определения антенны Ма для И =300 км равно

Рис• 1• Плотность распределения ошибки антенны P(8а) при h=300 км

Рис. 2. Плотность распределения ошибки антенны Р(да) при И=1000 км

9,6°, для И =1000 км равно 7,8°. Повышение точности с увеличением высоты объясняется уменьшением числа НС, затененных Землей.

Для удобства представления результатов была использована тройка углов ориентации (0у,р), задающая ориентацию СК

0Х1У222 относительно СК 0Х1Ї121. Система координат 0Х27222 может быть переведена в систему координат 0Х1У1тремя последовательными поворотами: 1) на угол у вокруг оси 0222; 2) на угол 0 вокруг оси 0У'2; 3) на угол р вокруг оси ОХ2, совпадающей с осью 0Х1.

Связь углов ориентации (0 у ,р) с найденным кватернионом V задается соотношениями [11]:

у 0

v0 = со8(—) cos(—) соз(-) + sin(—) sin(—) sin(—)

2 2

2

2 2

2

V1 = со'(у ) со'(0) ь'іп(р2) - ^п(у) ^п(у) ст(р2);

v2 = cos( у ) 8іп( 0 ) cos( Р ) + sin( у ) cos( 0 ) 8іп( Р );

22 у . . 0

22

2

v0 = sin(—)c0s(—)c0s(Р) - C0s(у )8іп(0)8іп(Р).

22

2

22

2

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

Изменение математических ожиданий

ошибок углов (у ,р) в зависимости от коэффициента о не превышает 0,3°.

Плотности распределения ошибок углов (0у,р), найденных по разработанному алгоритму, представлены на рис. 5-10.

Математические ожидания углов

(0 у ,р) при И=300 км

М[&,] = 2,8°, М{&¥]= 5,8°, М[5р] = 4,0;

при И =1000 км

М[50]= 1,5°, Ми ]= 5,1°, М5 ]= 3,6°.

О | 11 11111

^ -I-т-I-:—I—:-:—I—:-:---------------:-I-:-:-:-:-:-:-

12 4 6 8 11] 12 14 16 18 20

С

Рис. 3. Изменение математического ожидания ошибки 8в угла в от коэффициента а (И=300 км)

1---- -------------1-------------------^------^------;------£-----

12 4 6 3 10 12 14 16 18 20

С

Рис. 4. Изменение математического ожидания ошибки 8в угла в от коэффициента а (И=1000 км)

бе, градусы

Рис. 5. Плотность распределения Р (8в) ошибки угла в при И=300 км

и,и 0,11 0,1 0,09 0,08 0,07 0,06 0,05 0,04 0,03 0,02 0,01 0

5 градусы

Рис. 6. Плотность распределения Р (ду) ошибки угла у при И=300 км

40

0,26

0,24

0,22

0,2

0,18

0,16

0,14

0,12

0,1

0,08

0,06

0,04

0,02

О

-

V х

V X

V X N X \ X

V X

V X

V X \ X

V X

V X

V X

V X

вд

О 1

11111 КЗ СГ= 1

а= 10

I I I

9 10 11 12 13 14 15 16 17 18 19

5^, градусы

Рис. 7. Плотность распределения Р{89) ошибки угла (р при 11=300 км

Рис. 8. Плотность распределения Р (8в) ошибки угла в при И=1000 км

0,14 0,12 0,1 0,08 0,06 0,04 0,02 0

Рис. 9. Плотность распределения Р[§¥) ошибки угла у/ при И=1000 км

0 1

9 10 11 12 13 14 15 16 17 18 19

5 градусы

Рис. 10. Плотность распределения Р (5р) ошибки угла р при И=1000 км

Выводы

По результатам решения модельной задачи можно сделать следующие выводы.

1. Наибольший вклад в ошибку определения ориентации вносит ошибка определения вектора положения антенны в орбитальной системе координат. Снижение вклада данной ошибки возможно путем подбора коэффициента о.

2. С увеличением высоты полета погрешность определения ориентации КА уменьшается. Это объясняется уменьшением погрешности определения вектора положе-

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

Погрешность определения углов ориентации космического аппарата по предложенному алгоритму с вероятностью 90 % не превышает 5°.

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант РФФИ № 060-08-00244а).

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

1. Сидоров И. М., Прохоренко В. И. Определение углового положения искусственного спутника Земли с помощью датчиков магнитного поля // Космические исследования. 1968. - Т. VI. - Вып. 2. - С. 175-185.

2. Титов А. М., Антоненко В. В., Щукин В. П. Определение углового положения неориентированных ИСЗ по данным магнитометрических измерений // Космические исследования. - 1971. - Т. IX. - Вып. 3. -С.397-407.

3. Хацкевич И. Г. Определение ориентации ИСЗ по магнитометрическим измерениям //Космические исследования. - 1972. -Т. X. - Вып. 1. - С. 3-12.

4. Абрашкин В. И. и д.р. Определение вращательного движения спутника «Фотон-М2» по данным бортовых измерений магнитного поля Земли (Препринт Института прикладной математики им. М. в. Келдыша РАН, 2005, № 96).

5. Голубков В. В. Определение локальной ориентации космических аппаратов // Космические исследования. - 1970. - Т. VIII. -Вып. 6. - С. 811-822.

6. Титов А. М., Шукин В. П. Определение ориентации по двухвекторной системе измерений //Космические исследования. -1978. - Т. XVI. - Вып. 1. - С. 3-9.

7. Липтон А. Выставка инерциальных систем на подвижном основании. - М.: Наука, 1971.

8. Wertz J.R (Editor). Spacecraft Attitude Determination and Control. Dordrecht, The Netherlands. - 1978.

9. Wahba G. A Least Squares Estimate of Spacecraft Attitude //SIAM Review. - 1965., Vol.7, №3. - p. 409.

10. Коваленко А. П. Магнитные системы управления космическими летательными аппаратами. - М.: Машиностроение, 1976.

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

SPACE VEHICLE ATTITUDE CONTROL RECOVERY PROCEDURE COMBINING MAGNETOMETRIC AND PADIONAVIGATION MEASUREMENTS

© 2007 I. V. Belokonov, A. V. Kramlikh

Samara State Aerospace University

The paper discusses a procedure of space vehicle attitude control recovery combining magnetometric and radionavigation measurements. The efficiency of the procedure is confirmed on a model task.