CC BY
36
УДК 629.78 В01: 10.18287/2541-7533-2019-18-2-41-51
ОПРЕДЕЛЕНИЕ ДИНАМИКИ ВРАЩАТЕЛЬНОГО ДВИЖЕНИЯ КОСМИЧЕСКОГО АППАРАТА С ИСПОЛЬЗОВАНИЕМ ИНФОРМАЦИИ ГЛОБАЛЬНЫХ НАВИГАЦИОННЫХ СПУТНИКОВЫХ СИСТЕМ
© 2019
доктор технических наук, профессор, заведующий межвузовской кафедрой космических исследований; Самарский национальный исследовательский университет имени академика С.П. Королёва ibelokonov@mail.ru
кандидат технических наук, доцент, доцент межвузовской кафедры космических исследований; Самарский национальный исследовательский университет имени академика С.П. Королёва; kramlikh@mail.ru
аспирант;
Самарский национальный исследовательский университет имени академика С.П. Королёва igorlomaka63@gmail.com
Показана возможность применения навигационной аппаратуры потребителя в задаче определения динамики вращательного движения космического аппарата. Предложен подход, позволяющий оценивать параметры вращательного движения космического аппарата по анализу геометрической видимости навигационных космических аппаратов глобальных навигационных спутниковых систем. Предлагаемый подход заключается в аппроксимации накопленной информации о положении продольной оси космического аппарата с помощью модели углового движения и модели измерений. Аппроксимация производится исходя из минимизации суммы квадратов отклонений между вычисленными координатами вектора продольной оси и их смоделированными значениями. Процедура минимизации построена на алгоритме дифференциальной эволюции. Предложенный подход позволяет оценить угловые скорости космического аппарата с точностью не хуже 0.3 град/с и углы ориентации с точностью не хуже 15 градусов.
Космический аппарат; вращательное движение; навигация; дифференциальная эволюция.
Цитирование: Белоконов И.В., Крамлих А.В., Ломака И.А. Определение динамики вращательного движения космического аппарата с использованием информации глобальных навигационных спутниковых систем // Вестник Самарского университета. Аэрокосмическая техника, технологии и машиностроение. 2019. Т. 18, № 2. С. 41-51. DOI: 10.18287/2541-7533-2019-18-2-41-51
Введение
Для определения параметров вращательного движение (углов и угловых скоростей) космического аппарата (КА) традиционно используется информация от магнитометров [1], датчиков угловой скорости [2], звёздных датчиков [3], датчиков местной вертикали [4], датчиков направления на Солнце и т.д. Для решения задачи определения ориентации КА используется информация о токосъёме с панелей солнечных батарей [5]. Для обеспечения миссий КА, наряду с решением задачи определения параметров вращательного движения, требуется решить задачу навигации, т.е. определения параметров движения центра масс. Для решения задачи навигации на борту необходимо наличие навигационной аппаратуры потребителей (НАП), работающей по сигналам глобальных навигационных спутниковых систем (ГНСС) ГЛОНАСС (Россия), GPS (США), Galileo (Европейский союз). Современное состояние НАП позволяет устанавливать их на борт КА.
И. В. Белоконов
А. В. Крамлих
И. А. Ломака
Кроме решения задачи навигации НАП можно использовать для решения задачи определения ориентации. Так, например, НАП, использующая фазовые измерения, может быть применена для определения ориентации КА на основании известных принципов интерферометрии [6-8]. В настоящей работе предлагается подход, позволяющий использовать информацию, получаемую НАП, для оценки параметров вращательного движения КА.
Следует отметить, что не вся спутниковая радионавигационная информация в НАП используется в полном объёме. В работах [9;10] показано, что, используя информацию о положении навигационных космических аппаратов (НКА) ГНСС, о положении КА и геометрической видимости НКА со стороны КА, можно оценить пространственную ориентацию оси КА, на которой установлена антенна НАП. Погрешность оценки ориентации оси при этом не превысит 15° [9;10]. Достоинством данного подхода является возможность оценки пространственной ориентации оси КА по одномоментным измерениям. К недостаткам следует отнести невозможность получения оценки угловых скоростей КА, что является критическим для некоторых миссий.
Оценить параметры вращательного движения КА можно путём совместного использования информации о положении единичного вектора продольной оси КА и модели вращательного движения КА. Параметры вращательного движения КА (углы ориентации и угловые скорости) оцениваются при минимизации целевой функции, представляющей собой сумму квадратов разностей найденного единичного вектора продольной оси КА по алгоритму [9] и его модельных значений.
Пусть по результатам работы алгоритма [9] имеется набор измерений координат единичного вектора продольной оси КА в орбитальной системе координат (ОСК). Измерения получены на некотором интервале времени Т. Требуется получить оценки параметров углового движения КА на заданном интервале времени. Данная задача сводится к задаче нелинейной многопараметрической оптимизации, а именно к поиску минимума следующей целевой функции:
где (Ъ) - модель измерений координат единичного вектора продольной оси КА в ОСК; Ъ - вектор оцениваемых параметров; Т]сй - вычисленные координаты вектора
продольной оси.
Вектор оцениваемых параметров имеет вид
где - угол прецессии; а(^0) - угол атаки; <0) - угол собственного вращения;
а>х (¿{)),юу (¿0), сог (¿0) - угловые скорости, ¿0 - начальный момент времени.
Решение задачи проводится при следующих допущениях:
1) КА оснащён спутниковой навигационной антенной, обладающей диаграммой направленности в виде полусферы без обратных лепестков, что исключает возможность приёма сигнала с навигационных спутников, находящихся вне полусферы;
Постановка задачи
(1)
(2)
2) модуль угловой скорости КА не превышает 3 0 / с, что гарантирует получение полного кадра навигационной информации.
Математическая модель динамики движения КА
Для записи уравнений движения КА вокруг центра масс (точка О) (рис. 1), а также соотношений, используемых при обработке данных, вводятся три правые декартовые системы координат.
Рис. 1. Положение систем координат
Связанная с КА система координат (ССК) образована главными центральными осями инерции и имеет обозначение ОХсУс%с. Абсолютная геоцентрическая система
координат (АСК) имеет обозначение СХАУА2А с началом в центре масс Земли (точка С). Ось ХА направлена в точку весеннего равноденствия. Ось 2А направлена в северный полюс мира. Ось УА дополняет систему до правой. Орбитальная система координат (ОСК) имеет обозначение ОХ^0^0. Начало системы находится в центре масс КА. Ось Х0 направлена по радиус-вектору КА. Ось 20 перпендикулярна плоскости орбиты. Ось Уо дополняет систему до правой.
Переход от АСК к ОСК задаётся тремя последовательными поворотами на угол долготы восходящего узла ^ вокруг оси 2А, на угол наклонения орбиты I вокруг новой
оси ХА и на аргумент широты и вокруг новой оси 2А . Положение ССК относительно ОСК задаётся тремя последовательными поворотами на угол прецессии у вокруг оси У0, на угол атаки а вокруг новой оси 20 и на угол собственного вращения ф вокруг новой оси Уо .
Вращательное движение КА описывается динамическими уравнениями Эйлера и кинематическими уравнениями [11]:
Ях =М (yVz - Va2Xa3l ) + pE (ryVz - ГУсу ) S
1 — X X
Vy =-(x®z — Va2ia31 ) + ^777 pE (х — rxVcz ) S I Vc I ,
1Л \ Л Z L1 jl /л >1
+ X/ 1 + X/
d> z = — (1 — X + X/(w — va2ia31 ) + XpE (ryVcz — rVy ) S
(3)
4 = 05 (—(х —®г)Ч1 —(®y —Vr ) 42 —(z —Vr )Чз ),
4 = 0.5 (( —+(z —Vr И —(y —Vr )Чз ),
q = °.5 (( — G)r)qo +(x —ar )чз —(o)2 —ar )q ), 4 = 0.5 (( — ®rK +(y —Vr )q —(х —Vr )42).
(4)
Здесь Х = 1х/1г и ¡л = (1у -х - безразмерные коэффициенты инерции (выражения для безразмерных коэффициентов инерции взяты из [1]); у = 3^,е/г 3 - коэффициент гравитационного момента; р = 0.5рСх - коэффициент аэродинамического момента; S -площадь проекции поверхности КА на плоскость, перпендикулярную вектору скорости; р - плотность атмосферы на высоте орбиты КА; Ус = (УСх ¥у УС2) - вектор орбитальной скорости в ССК; Г =( гх гу т2) - вектор, направленный из центра масс в
центр давления в ССК; Сг = Л(д) |0 0 сорб ] сорб - угловая орбитальная скорость КА.
Матрица А перехода из ОСК в ССК имеет следующий вид:
A(q) =
1 — 2 (42 + Чз Чз) 2 (4142 + ЧзЧо) 2 ( 41 Чз + 42 Чо) 2 ( + Чз Чо) 1 — 2 (Ч4 + ЧзЧз) 2 ( Чз + 44о)
2 (qq + 42 Чо) 2 ( Чз + 41 Чо) 1—2 (qq + 42 42)
Модель измерений
Алгоритм определения ориентации оси КА основывается на использовании информации о пространственном положении НКА ГНСС ГЛОНАСС и GPS [9]. Для определённости будем считать, что антенна расположена по продольной оси КА. Задача определения ориентации продольной оси КА сводится к отысканию оценки вектора
направляющих косинусов фазового центра антенны НАП в ОСК Ao =(хо yo zo )Г , расположенной по продольной оси КА, из условия минимума целевой функции Ф (хо Уо zo), отражающей условия видимости/невидимости НКА:
cos (ao, gradBi )> cos (a), (i = 1, NB),
cos (ao, gradHBJ )> cos (a), (J =1 NHB ),
где gradi =( xoi, yoi, zoi) - единичный вектор дальности до z'-го НКА в проекциях на оси ОСК; NB, NHB — количество видимых и невидимых НКА соответственно; а - угол
полураствора конуса затенения.
Процедура решения задачи определения ориентации продольной оси КА включает следующие этапы:
1. Расчёт эфемерид невидимых навигационных спутников ГНСС ГЛОНАСС и GPS на моменты времени решения задачи определения ориентации.
2. Пересчёт дальностей до видимых/невидимых НКА из АСК в ОСК.
3. Исключение из рассмотрения невидимых НКА, затенённых Землёй. Условия затенения Землёй имеют вид:
Г
Z2k < 0 и Ы > C0S
arcsin
f R ЛЛ
v r + h уу
k =
1 N
I 5 глонасс5
1 N
l ' GPS'
4. Отыскание оценки вектора направляющих косинусов фазового центра антенны из условия минимума целевой функции, отражающей условия видимости/невидимости НКА (рис. 2):
Ф (х у г ) = "V (х х + у у + г г -1) (х х + у у + г г -1) , (5)
V о У о о } / ' \ ог о У ог У о ог о / / < \ о] о У о] У о о] о у' V/
г=1 г=1
где ЫВ - количество видимых НКА; ЫНВ - количество невидимых НКА.
Кшдш ииЛ я екгс*. НК i-ЦА
ItiTHHHUi НК i ДА
Едиккчниг .илпряикиил ИЛ ИЩИ ta гг НКА
Пл f-; к-) с :v эятгксшгл 1«кгт|>ухикеК КА
Нд 11 нич и Мс и ¡шрам! нкя Н1НСЕ11ЛИКЫГ Н VA
Рис. 2. Иллюстрация геометрической видимости/невидимости НКА
Процедура минимизации целевой функции (5) сводится к решению системы трёх линейных алгебраических уравнений:
f N, N„, Л f N, Nя. Л f N к, Л N, K,
Zx° + Z ХоУ xo + Z x y ■ +Z x y ■ / , oiijoj / , otj oi yo + Z xoizoj + Z xoizoj zo = Z x2 + Z x=,
V i = j= V i=1 j=1 У V i =1 j=1 У i=1 j=1
f Ne NH„ Л f N, NH, Л f N Л N, Nн,
Z x y . + Z x / ' ois oj / ' ъУ oj xo + Z Z У oi ^ Z 1 У o2 yo + Z У -z ■ / * ^ oi oj + Z yoizoj zo Z \ У oi +Z у2.
v ч i=1 j=1 v1=1 j=1 У V i=1 j=1 У i=1 j=1
f N, N Л f N Л f N Л N, Nн,
Z x .z ■ + Z x .z ■ / ' oi oi / , oi oj xo + Z У -z ■ +Z У -z ■ / * ^ oi oj ¿^ J oi oj yo + Z 2 , 2 zoi +Z zoj zo = Z zo2i + Z zj.
V i =1 j=1 V i=1 j=1 У V i=1 j=1 У i=1 j=1
Использование данного алгоритма позволяет получить ориентацию оси КА с точностью не хуже 3а = 15 град.
Подход к определению параметров углового движения КА
Предлагаемый подход заключается в минимизации целевой функции (1) с использованием моделей движения (3), (4) и моделей измерений (5), (6). В результате минимизации (1) определяется вектор оцениваемых параметров (2). Подход состоит из следующих этапов:
1) в течение промежутка времени T с помощью соотношений (5), (6) вычисляются координаты вектора продольной оси КА. Полученные координаты накапливаются и образуют массив
2) полученный массив измерений обрабатывается с помощью алгоритма дифференциальной эволюции (АДЭ), в результате чего происходит минимизация функции (1) и определяется вектор оцениваемых параметров (2).
Численный метод решения задачи
Процедура поиска минимума функции (1) построена на применении АДЭ [12], который заключается в выполнении следующих этапов:
1) задание исходного массива X, состоящего из N векторов Ь (N на каждой итерации постоянно и является одним из параметров АДЭ), причём каждый элемент вектора Ъ генерируется случайным образом в требуемых пределах;
2) генерация нового массива Хтк векторов Ь : для каждого вектора Ь из массива
X выбираются случайным образом три различных вектора Ь1з Ь2, Ь3 , индексы которых не совпадают с индексом вектора Ь , и вычисляется так называемый модифицированный вектор: Ьm = Ь1 + F(Ь2 - Ь3), где F - константа в интервале [0 1], являющаяся параметром алгоритма АДЭ;
3) над модифицированным вектором Ьп1 выполняется операция «скрещивания», состоящая в том, что некоторые его элементы замещаются соответствующими элементами из исходного вектора Ь (каждый элемент замещается с вероятностью P, которая также является параметром алгоритма);
4) если полученный вектор оказывается лучше вектора Ь, то есть значение целевой функции уменьшается: J(Ь{ ) < J(Ьi), то в новом массиве Хпек вектор Ь заменяется на новый вектор Ь( (пробный вектор), в противном случае сохраняется Ь;
5) пункты 1 - 4 повторяются до тех пор, пока не будет выполнен критерий остановки.
Критерий остановки зависит от решаемой задачи. В данной работе итерационная процедура останавливается, как только разность между значениями целевой функции на текущей и предыдущей итерациях становится меньше 0,001.
Статистическое исследование точности решения задачи восстановления
вращательного движения КА
Для оценки точности определения углового движения было проведено статистическое моделирование решения задачи.
Моделирование происходило следующим образом:
1. Было сгенерировано 400 вариантов начальных условий углового движения ^х (^),®у (^)(^)),«(^),ф(t0) в предположении равновероятного закона их
распределения. Полагалось, что > 3 град/с, ^(?0) и ^(?0) лежат в диапазоне от 0
до 360 градусов, а(?0) лежит в диапазоне от 0 до 180 градусов.
2. Для каждого варианта начальных условий были получены координаты вектора продольной оси КА с точностью не хуже 3а = 15 град на интервале времени T = 1500 с.
3. Результаты решения задачи сравнивались с начальными условиями, которые использовались для моделирования измерений. В результате были получены ошибки алгоритма восстановления вращательного движения.
Результаты моделирования приведены на рис. 3, 4, где использованы следующие обозначения: М - математическое ожидание, а - среднеквадратическое отклонение.
Рис. 3. Распределение вероятности ошибки по угловой скорости
-5 0 5
Ошибку, град
-5 0 5
Ошибка а, град
0 -20
-5 0
Ошибка
град
Рис. 4. Распределение вероятности ошибки по углам ориентации
Заключение
Предложенный подход позволяет оценить три угла ориентации и три угловые скорости КА.
Наихудшая точность в определении ориентации составила а — 5 град, а по угловой скорости составила а — 0,12 град/с.
Исследование выполнено за счёт гранта Российского научного фонда (проект №17-79-20215).
Библиографический список
1. Абрашкин В.И., Воронов К.Е., Пияков А.В., Пузин Ю.Я., Сазонов В.В., Сем-кин Н.Д., Филиппов А.С., Чебуков С.Ю. Определение вращательного движения спутника АИСТ по данным бортовых измерений магнитного поля Земли // Препринты ИПМ им. М.В. Келдыша. 2014. № 17. 38 с.
2. Беляев М.Ю., Волков О.Н., Монахов М.И., Сазонов В.В. Оценка точности методики реконструкции вращательного движения спутника по измерениям его угловой скорости и магнитного поля Земли // Препринты ИПМ им. М.В. Келдыша. 2016. № 4. 40 с. DOI: 10.20948/prepr-2016-4
3. Аванесов Г.А., Бессонов Р.В., Куркина А.Н., Никитин А.В., Сазонов В.В. Определение движения космического аппарата по измерениям четырёх звёздных датчиков // Препринты ИПМ им. М.В. Келдыша. 2016. № 57. 38 c. DOI: 10.20948/prepr-2016-57
4. Lomaka I., Belokonov I., Ustugov E. Technology for determining the local vertical of nanosatellite by processing videoimages of the Earth horizon // IFAC-PapersOnline. 2016. V. 49, Iss. 17. P. 206-211. DOI: 10.1016/j.ifacol.2016.09.036
5. Беляев М.Ю., Матвеева Т В., Монахов М.И., Сазонов В В., Цветков В.В. Определение вращательного движения кораблей Прогресс по данным измерений угловой скорости и токосъёма с солнечных батарей // Препринты ИПМ им. М.В. Келдыша. 2012. № 39. 36 c.
6. Степанов О.А, Кошаев Д.А. Исследование методов решения задачи ориентации с использованием спутниковых систем // Гироскопия и навигация. 1999. № 2 (25). C. 30-55.
7. Treder A. Attitude error from phase center uncertainty in an interferometric GPS antenna array // Collection of Technical Papers - AIAA/AAS Astrodynamics Specialist Conference (August, 21-24, 2006, Keystone, CO; United States). 2006. V. 2. P. 904-917. DOI: 10.2514/6.2006-6384
8. Torre A.D., Caporali A., Praticelli N., Facchinetti C. Attitude and direction sensor using GPS carrier phase data // European Space Agency, (Special Publication) ESA SP. 2004. V. 548. P. 291-296.
9. Белоконов И.В., Крамлих А.В. Методика восстановления ориентации космического аппарата при комплексировании магнитометрических и радионавигационных измерений // Вестник Самарского государственного аэрокосмического университета имени академика С П. Королёва. 2007. № 1 (12). C. 22-30.
10. Григорьева М.Е., Крамлих А.В. Адаптивный алгоритм определения ориентации низковысотных космических аппаратов на основе обработки одномоментных разнотипных измерений // Вестник Самарского государственного аэрокосмического университета имени академика С.П. Королёва (национального исследовательского университета). 2012. № 4 (35). C. 69-75.
11. Белецкий В.В. Движение искусственного спутника относительно центра масс. М.: Наука, 1965. 416 с.
12. Storn R., Price K. Differential evolution - a simple and efficient heuristic for global optimization over continuous spaces // Journal of Global Optimization. 1997. V. 11. P. 341359. DOI: 10.1023/A:1008202821328
DETERMINING THE DYNAMICS OF SPACECRAFT ROTATIONAL MOTION WITH THE USE OF INFORMATION PROVIDED BY GLOBAL NAVIGATION SATELLITE SYSTEMS
© 2019
I. V. Belokonov Doctor of Science (Engineering), Professor, Head of Inter-University Department of Space Research;
Samara National Research University, Samara, Russian Federation; ibelokonov@mail.ru
A. V. Kramlikh Candidate of Science (Engineering), Associate Professor, Inter-University Department of Space Research;
Samara National Research University, Samara, Russian Federation; kramlikh@mail.ru
I. A. Lomaka Postgraduate Student;
Samara National Research University, Samara, Russian Federation; igorlomaka63@gmail.com
The paper shows the possibility of using navigation user equipment in the task of determining the dynamics of spacecraft rotational motion. An approach is proposed that makes it possible to estimate the parameters of the rotational motion of a spacecraft by analyzing the geometric visibility of navigational spacecraft of global navigation satellite systems. The proposed approach consists in approximation of the accumulated information on the position of the spacecraft longitudinal axis using the angular motion model and the measurement model. The approximation is made on the basis of minimizing the sum of squared deviations between the calculated coordinates of the vector of the longitudinal axis and their simulated values. The minimization procedure is based on the algorithm of differential evolution. The proposed approach allows us to estimate the angular velocity of the nanosatellite with an accuracy of at least 0.3 deg/s and orientation angles with an accuracy of at least 15 degrees.
Spacecraft; rotational motion; navigation; differential evolution.
Citation: Belokonov I.V., Kramlikh A.V., Lomaka I.A. Determining the dynamics of spacecraft rotational motion with the use of information provided by global navigation satellite systems. Vestnik of Samara University. Aerospace and Mechanical Engineering. 2019. V. 18, no. 2. P. 41-51. DOI: 10.18287/2541-7533-2019-18-2-41-51
References
1. Abrashkin V.I., Voronov K.E., Piyakov A.V., Puzin Yu.Ya., Sazonov V.V., Sem-kin N.D., Filippov A.S., Chebukov S.Y. Determination of the spacecraft AIST attitude motion on measurements of the Earth magnetic field. Keldysh Institute Preprints. 2014. No. 17. 38 p. (In Russ.)
2. Belyaev M.Yu., Volkov O.N., Monakhov M.I., Sazonov V.V. Accuracy estimation of the technique for reconstructing the spacecraft attitude motion by measurements of its angular rate and Earth magnetic field strength. Keldysh Institute Preprints. 2016. No. 4. 40 p. DOI: 10.20948/prepr-2016-4. (In Russ.)
3. Avanesov G.A., Bessonov R.V., Kurkina A.N., Nikitin A.V., Sazonov V.V. Determination of a spacecraft attitude motion by measurements of four star sensors. Keldysh Institute Preprints. 2016. No. 57. 38 p. DOI: 10.20948/prepr-2016-57. (In Russ.)
4. Lomaka I., Belokonov I., Ustugov E. Technology for determining the local vertical of nanosatellite by processing videoimages of the Earth horizon. IFAC-PapersOnline. 2016. V. 49, Iss. 17. P. 206-211. DOI: 10.1016/j.ifacol.2016.09.036
5. Belyaev M.Yu., Matveeva T.V., Monakhov M.I., Sazonov V.V., Tsvetkov V.V. Reconstruction of spacecraft Progress attitude motion by measurements of their angular rate and electric current from solar batteries. Keldysh Institute Preprints. 2012. No. 39. 36 p. (In Russ.)
6. Stepanov O.A., Koshaev D.A. Analysis of solution methods for orientation problem
using satellite systems. Giroskopiya i Navigatsiya. 1999. No. 2 (25). P. 30-55. (In Russ.)
7. Treder A. Attitude error from phase center uncertainty in an interferometric GPS antenna array. Collection of Technical Papers - AIAA/AAS Astrodynamics Specialist Conference (August, 21-24, 2006, Keystone, CO; United States). 2006. V. 2. P. 904-917. DOI: 10.2514/6.2006-6384
8. Torre A.D., Caporali A., Praticelli N., Facchinetti C. Attitude and direction sensor using GPS carrier phase data. European Space Agency, (Special Publication) ESA SP. 2004. V. 548. P. 291-296.
9. Belokonov I.V., Kramlikh A.V. Space vehicle attitude control recovery procedure combining magnetometric and padionavigation measurements. Vestnik of the Samara State Aerospace University. 2007. No. 1 (12). P. 22-30. (In Russ.)
10. Grigoryeva M.Ye., Kramlikh A.V. Adaptive algorithm of determining low altitude spacecraft orientation on the basis of processing instant diverse-type measurements. Vestnik of the Samara State Aerospace University. 2012. No. 4 (35). P. 69-75. (In Russ.)
11. Beletskiy V.V. Dvizhenie iskusstvennogo sputnika otnositel'no tsentra mass [Motion of an artificial satellite relative to the center of mass]. Moscow: Nauka Publ., 1965. 416 p.
12. Storn R., Price K. Differential evolution - a simple and efficient heuristic for global optimization over continuous spaces. Journal of Global Optimization. 1997. V. 11. P. 341359. DOI: 10.1023/A:1008202821328
