«Баллистико-навигационное проектирование полётов к Луне, планетам и малым телам Солнечной системы ...»

УЧРЕЖДЕНИЕ РОССИЙСКОЙ АКАДЕМИИ НАУК

Институт прикладной математики им. М.В. Келдыша РАН

На правах рукописи

ТУЧИН Андрей Георгиевич

Баллистико-навигационное

проектирование полётов к Луне, планетам

и малым телам Солнечной системы

Специальность 01.02.01 – Теоретическая механика

Диссертация

на соискание учёной степени доктора физико-математических наук

Москва – 2010 Содержание Обозначения и сокращения

Введение

Глава 1 Проектирование квазисинхронных орбит КА вокруг Фобоса для решения задачи посадки на его поверхность... 18 1.1 Численный алгоритм построения множества КСО с минимальным дрейфом.. Формализация постановки задачи

1.1. Упрощённая модель движения КА относительно Фобоса

1.1. Область перебора

1.1. Переменные перебора

1.1. Численный анализ

1.1. Результаты численного анализа

1.1. Сравнение характеристик движения в полной и упрощённой моделях.................. 1.1. 1.2 Начальное приближение для расчёта множества КСО

1.2.1 Упрощённые уравнения движения

1.2.2 Предварительный численный анализ. Постановка задачи аналитических исследований

1.2.3 Подготовка системы (1.19-1), (1.19-2) к усреднению

1.2.4 Частный случай плоского кругового движения задачи Хилла ( = 1)

d

1.2.5 Вычисление значения интеграла r I I I

1.2.6 Вычисление производных,, A a 1.2.7 Усреднённые уравнения движения

1.2.8 Условие для поиска начального приближения начальных условий.

Выбор параметра

1.2.9 Алгоритм формирования начального приближения

1.2.10 Вычислительная процедура поиска минимума

1.2.11 Результаты расчётов

1.2.12 Алгоритм проектирования КСО

Глава 2 Баллистика, навигация и управление движением КА на этапе его посадки на поверхность Фобоса

2.1 Сближение с Фобосом и посадка на его поверхность: общая схема

2.1.1 Сближение с Фобосом

2.1.2 Условия посадки

2.1.3 Навигационные приборы обеспечения посадки

2.1.4 Схема управляемой посадки

2.1.5 Управление включением двигателей

2.2 Алгоритмы управления движением КА

2.2.1 Бортовые алгоритмы навигации и управления

2.2.2 Уравнения движения КА относительно Фобоса

2.2.3 Расчёт коррекций на участке перелёта от момента схода с КСО до точки начала вертикального спуска

2.2.4 Определение вектора состояния КА по измерениям лазерного высотомера и доплеровской системы

2.3 Реализация алгоритмов посадки в среде операционной системы реального времени

Глава 3 Определение параметров движения КА по результатам измерений при наличии немоделируемых ускорений.......... 3.1 Алгоритм оценки вектора состояния и суммарных воздействий возмущений между измерениями

3.1.1 Постановка задачи

3.1.2 Линейный случай

3.1.3 Нелинейный случай

3.1.4 Проверка качества измерений с использованием приведённого среднеквадратического отклонения

3.2 Алгоритм оценки вектора состояния в случае отсутствия шума

3.3 Алгоритм оценки вектора состояния и средних значений приращений возмущений

3.4 Алгоритм оценки вектора состояния с использованием метода мешающих параметров

3.4.1 Мешающие параметры в форме белого шума

3.4.2 Мешающие параметры в форме случайных величин, постоянных на всем интервале

3.5 Оценка вектора состояния и возмущений дискретной динамической системы и свойства этих оценок

Глава 4 Вопросы баллистики и навигации в проектах полётов к Луне, планетам и малым телам Солнечной системы..... 4.1 Проблемы баллистики и навигации в проектах полётов к Луне, точке L2 системы Солнце – Земля, планетам и малым телам Солнечной системы

4.1.1 Проект «Фобос-Грунт»

4.1.2 Российские проекты полёта к Луне

4.1.3 Полёт к точке L2 системы Солнце – Земля

4.1.4 Проект полёта к Венере с целью доставки на её поверхность долгоживущей станции, баллонов в атмосферу Венеры и выхода орбитальной станции на орбиту искусственного спутника Венеры

4.1.5 Проект полёта к Юпитеру с целью посадки КА на поверхность Европы, естественного спутника Юпитера (проект «Лаплас»)

4.1.6 Особенности баллистико-навигационного обеспечения полётов к дальним планетам

4.1.7 Подготовка и проведение гравитационного маневра КА «Розетта»

4.1.8 Проект полёта к Солнцу

4.2 Методы обеспечения навигации и управления КА дальнего космоса

4.2.1 Управление и навигация на участках перелёта к планетам с использованием ЖРД

4.2.2 Управление и навигация на участках перелёта к планетам с использованием ЭРДУ

4.2.3 Управление и навигация при проведении гравитационных маневров

4.2.4 Автономная навигация по изображениям планет на участках перелёта от Земли к планете-цели

4.2.5 Управление и навигация на орбитах искусственных спутников планет

4.3 Обобщённая баллистическая задача

Глава 5 Анализ точности приведения КА к Луне и планетам........ 5.1 Методика и вычислительный алгоритм оценки точности определения и прогнозирования параметров движения КА по наземным траекторным измерениям

5.1.1 Расчёт ковариационной матрицы ошибок определения вектора состояния, обусловленных неучтёнными ускорениями и случайными составляющими ошибок измерений

5.1.2 Поправки на ошибки, обусловленные систематическими составляющими ошибок измерений

5.1.3 Поправки на ошибки, обусловленные ошибками привязки фазового центра антенны

5.1.4 Поправки на ошибки, обусловленные тропосферными составляющими ошибок измерений

5.1.5 Поправки на ошибки, обусловленные ионосферными составляющими ошибок измерений

5.1.6 Прогнозирование ковариационной матрицы

5.2 Методика выбора схемы проведения коррекций при перелётах к Луне и планетам

5.3 Примеры расчётов определения точности приведения КА к Луне и планетам с использованием жидкостных реактивных двигателей................ 5.3.1 Определение точности приведения КА к Марсу (на примере проекта «Фобос-Грунт»)

5.3.2 Определение точности приведения КА к Луне (на примере проекта «Луна-Глоб»)

5.3.3 Определение точности приведения КА к Венере (на примере проекта «Венера-Д»)

5.4 Определение точности приведения КА с ЭРДУ к Юпитеру в проекте «Лаплас»

5.4.1 Моделирование определения орбиты

5.4.2 Модель ошибок ЭРДУ

5.4.3 Расширенная система

5.4.4 Оценка ошибок определения вектора состояния и параметров линейных функций

5.4.5 Технология управления КА на участке перелёта

5.4.6 Расчёт ошибок приведения КА к Юпитеру и оценки дополнительного расхода топлива

5.5 Требования по точности наземных и бортовых траекторных измерений........ 5.5.1 Требования к точности наземных траекторных измерений

5.5.2 Требования к точности бортовых акселерометров

5.5.3 Требования к точности бортовых измерений направления на планеты, их спутники и астероиды

Заключение

Список использованных источников

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

Система обозначений Для скалярных переменных используется курсивный шрифт, для векторов и матриц – прямой шрифт, при этом матрицы обозначены прописными буквами, а векторы – строчными. Исключение составляют греческие буквы, которые всегда представлены в прямом начертании. Матрица для преобразования из системы координат p в систему координат q обозначается, как Cq.

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

Перечень сокращений CCSDS – The Consultative Comitee for Space Data Systems;

Delta-DOR – Delta Differencial One-way Ranging measurement;

DSN – Deep Space Network;

ESA – European Space Agency;

GPS – Global Positioning System;

RNB – Radial – Normal (=Transversal) – Binormal reference frame;

АОНС – автономная оптическая навигационная система;

АСП – автономная система посадки;

БИБ – бесплатформенный инерциальный блок;

БКУ – бортовой комплекс управления;

БНО – баллистико-навигационное обеспечение;

БОКЗ-МФ – блок определения координат звёздный, модифицированный;

ВА – возвращаемый аппарат;

ДИСД – доплеровский измеритель скорости и дальности;

ДМТ – двигатель малой тяги;

ДУ – двигательная установка;

ЖРД – жидкостный реактивный двигатель;

ЗГ – задающий генератор;

ИП – измерительный пункт;

ИСЗ – искусственный спутник Земли;

ИСМ – искусственный спутник Марса;

КА – космический аппарат;

КСО – квазисинхронная орбита;

ЛВ – лазерный высотомер;

ЛКИ – лётно-конструкторские испытания;

МДУ – маршевая двигательная установка;

НТО – научно-технический отчёт;

ОИСМ – орбита искусственного спутника Марса;

ОНА – остронаправленная антенна;

ПА – посадочный аппарат;

ПМ – перелётный модуль;

ПСК – приборная система координат;

ПСО – постоянная солнечная ориентация;

РН – ракета-носитель;

РСДБ – радиоинтерферометрические измерения со сверхдлинной базой;

РЭ – референц-эллипсоид;

СК – система координат;

СКО – среднеквадратичная ошибка;

ТМИ – телеметрическая информация;

ЦПФ – цифровая поверхность Фобоса;

ЭМИО – электромаховичные исполнительные органы;

ЭРДУ – электроракетная двигательная установка.

Сокращения в индексах E – Земля (Earth);

grad – градиент;

obs – измеренный (observed) Ph – Фобос (Phobos);

prog – программный;

rcv – приёмник (receiver);

real – реальный;

RF – вращающаяся система координат (rotating frame);

S – Солнце (Sun);

sc – космический аппарат (spacecraft);

snd – передатчик (sender);

td – «сухая» (dry) поправка;

tw – «мокрая» (wet) поправка.

Введение Диссертация посвящена теоретико-механическим вопросам проектирования полетов космических аппаратов (КА) к планетам Солнечной системы и их естественным спутникам. Разработан метод проектирования квазисинхронных орбит КА вокруг Фобоса, предназначенных для обеспечения посадки КА на его поверхность. Разработаны методы определения параметров движения КА по траекторным измерениям на фоне работы двигательной установки; разработаны методы расчета точности приведения КА к планете цели при использовании химических и электроракетных двигателей с учетом ошибок прогноза движения КА, ошибок исполнения маневров и коррекций и различных возмущений, вносимых в движение центра масс КА. Разработанные методы опираются на опыт баллистиконавигационного обеспечения полетов к Венере КА «Венера-15,16», «Вега-1,2», к Марсу КА «Фобос-2» и нашли применение в проекте «Фобос-Грунт», при проектировании полетов к Луне в проектах «Луна-Ресурс» и «Луна-Глоб», к Венере (проект «Венера-Д»), в систему Юпитера с целью посадки КА на его естественный спутник Европу. Это определяет актуальность и практическую значимость диссертации.

Цель работы состоит в разработке теоретико-механических и математических методов, обеспечивающих баллистику и навигацию в проекте «Фобос-Грунт» и распространению этих методов на решение задач навигации и управления в перспективных отечественных проектах полётов к Луне (проекты «Луна-Ресурс» и «Луна-Глоб»), к Венере (проект «Венера-Д»), и в систему Юпитера (проект «Лаплас»).

Результаты, представленные в диссертации, докладывались автором на 16-м симпозиуме IFAC по автоматическому управлению (Санкт-Петербург, июнь 2004 г);

на 17-м (Москва, июнь 2003г) и 18-м (Мюнхен, октябрь 2004г.) Международных симпозиумах по динамике космического полёта, на Общероссийском семинаре «Современные методы навигации и управления движением» (Институт проблем Управления РАН, 31 марта 2009 г.), на семинаре «Солнечная система и смежные проблемы физики и механики» (ИПМ им. М.В. Келдыша РАН, 18 марта 2008 г.).

Список основных публикаций по теме диссертации состоит из 18 работ; работ из этого списка опубликованы в изданиях, рекомендованных ВАК; 5 из них в соавторстве. Из работ, выполненных с соавторами, в диссертацию включены только результаты, полученные автором.

Основные результаты, полученные в диссертации, приведены в заключении.

В первой главе рассмотрены методы проектирования квазисинхронных орбит КА вокруг Фобоса для решения задачи посадки на его поверхность. Результаты, представленные в этой главе, описаны в работах: [112 – Tuchin, 2003; 57 – Тучин, 2007; 60 – Тучин, 2008; 61 – Тучин, 2009; 59 – Тучин, 2008] непосредственного сближения КА с Фобосом. На этом этапе должны быть получены навигационные данные, необходимые для успешной работы автономной системы посадки. Для решения этой задачи необходима такая орбита КА, двигаясь по которой КА не удалялся бы более, чем на 70 км от Фобоса. К классу таких орбит в терминологии ограниченной задачи трёх тел относятся квазиспутниковые орбиты [37 – Лидов, 1994]. Это орбиты с обратным движением, охватывающие тело меньшей массы и расположенные вне сферы его действия. Ниже будем называть такие орбиты квазисинхронными. Двигаясь по такой орбите, КА облетает Фобос за интервал времени, который примерно равен периоду обращения Фобоса вокруг Марса.

Моделирование показало, что требуются четыре дня полёта КА на КСО для определения параметров движения с необходимой точностью [67 – Тучин, 2002;

66 – Тучин, 2010]. Поэтому интервал пребывания на КСО оценивается пятью сутками. К номинальной КСО предъявляется ряд требований. За время пребывания на КСО с борта КА должно быть получено телевизионное изображение предполагаемого района посадки и проверена работоспособность лазерного высотомера автономной системы посадки. В момент начала сеанса посадки, при получении изображения района посадки и проверки работоспособности лазерного высотомера, удаление от поверхности Фобоса не должно превосходить 60 км и должны выполняться условия по освещённости и радиосвязи со станциями слежения в Евпатории и Уссурийске. Если в штатном сеансе работы автономной системы посадка не будет осуществлена и КА останется на прежней орбите, должна быть обеспечена возможность повторения попытки посадки на Фобос.

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

Вторая глава посвящена алгоритмам навигации и управления в схеме посадки на Фобос. Результаты, представленные в этой главе, описаны в работах [114 – Tuchin, 2004; 63 – Тучин, 2009]. В этой главе содержится модель движения КА относительно Фобоса, которая достаточно проста для бортовых вычислений, обеспечивает требуемую точность и позволяет рассчитывать коррекции и решать навигационные задачи. Далее строится алгоритм оценки вектора состояния КА по измерениям лазерного высотомера и доплеровского измерителя скорости и дальности с использованием цифровой модели поверхности Фобоса, разработанной специалистами Института геохимии и аналитической химии им. В.И. Вернадского РАН. В этой главе анализируется схема посадки и условия ее выполнения.

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

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

определение параметров движения КА с электроракетной двигательной установкой (ЭРДУ).

К вопросам, которые рассматриваются в диссертации, в первую очередь, относятся проблемы решения задач навигации и управления полётом с включённым ЭРДУ. Фактическое ускорение, создаваемое ЭРДУ, отличается от модели этого ускорения, заложенного в расчёты. Имеются ошибки величины и ориентации вектора тяги ЭРДУ в пространстве.

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

Результаты, представленные в данной главе, описаны в работе [62 – Тучин, 2009] в части алгоритмов определения параметров движения КА и в работе [58 – Тучин, 2004].

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

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

метод наименьших квадратов;

расширенный фильтр Калмана.

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

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

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

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

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

Минимизация функционала выполняется итерационно. На каждом шаге итерационного процесса определяется поправка к искомым параметрам.

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

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

Четвертая глава содержит анализ проблем баллистики и навигации в проектах полётов к Луне, точке L2 системы Солнце – Земля и планетам Солнечной системы.

Рассмотрены отечественные проекты «Фобос-Грунт», полёты к Луне, в окрестность точки L2 системы Солнце – Земля, к Венере и Юпитеру, а также проекты «Кассини», «New Horizons» и «Пионер-10, -11» (США).

соответствуют следующей обобщённой задаче.

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

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

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

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

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

Предлагаемый метод гарантированной оценки точности решения обобщённой баллистической задачи рассмотрен в пятой главе.

Результаты, представленные в пятой главе, описаны в разделах анализа точности определения и прогнозирования параметров движения КА в проекте «Фобос-Грунт» [61 – Тучин, 2009; 66 – Тучин, 2010; 67 – Тучин, 2002, 113 – Tuchin, 2004]. Эти методы применимы для анализа точности определения и прогнозирования параметров движения КА в полётах к Луне, планетам и малым телам Солнечной системы и были использованы при подготовке научнотехнических отчётов [122 – НТО 5-14-04; 124 – НТО 5-04-06; 132 – НТО 5-31-07] по проекту «Фобос-Грунт»; [128 – НТО 5-06-08; 135 – НТО 5-02-09] по проектам полёта к Луне; [141 – НТО 5-013-09] по проекту «Венера-Д».

Фактическая точность приведения КА к планете-цели определяется ошибками прогнозирования параметров движения, ошибками исполнения маневров или коррекций, а также неучтёнными возмущениями в движении центра масс КА.

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

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

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

Выполняя расчёты по описанному выше методу, получаем характеристики рассматриваемой схемы проведения коррекций: точность приведения и затраты характеристической скорости. Далее, варьируя числом коррекций и моментами их исполнения, получаем приемлемую схему.

Для применения метода необходимо иметь функциональную связь ошибок по целевым параметрам с отклонением вектора скорости от его значения по опорной траектории. Например, при решении задачи приведения КА к Луне или планете – эта функциональная связь определяется зависимостью значений параметров картинной плоскости от компонент вектора скорости. При решении задач приведения КА на гало-орбиту вокруг точки L2 эта зависимость определяется зависимостью от вектора скорости параметров, определяющих геометрию галоорбиты.

Работа метода показана на задачах приведения КА к Луне или планете.

Глава 1 Проектирование квазисинхронных орбит КА вокруг Фобоса для решения задачи посадки на его поверхность Эта глава посвящена вопросам выбора множества квазисинхронных орбит (КСО), с которых обеспечивается посадка на Фобос. Результаты, представленные в этой главе, описаны в работах [112 – Tuchin, 2003; 57 – Тучин, 2007; 60 – Тучин, 2008; 61 – Тучин, 2009; 59 – Тучин, 2008].

непосредственного сближения КА с Фобосом. На этом этапе должны быть получены навигационные данные, необходимые для успешной работы автономной системы посадки. Для решения этой задачи необходима такая орбита КА, двигаясь по которой КА не удалялся бы более, чем на 70 км от Фобоса. К классу таких орбит в терминологии ограниченной задачи трёх тел относятся квазиспутниковые орбиты [37 – Лидов, 1994]. Это орбиты с обратным движением, охватывающие тело меньшей массы и расположенные вне сферы его действия. Ниже будем называть такие орбиты квазисинхронными. Двигаясь по такой орбите, КА облетает Фобос за интервал времени, который примерно равен периоду обращения Фобоса вокруг Марса.

Моделирование показало, что требуются четыре дня полёта КА на КСО для определения параметров движения с необходимой точностью [67 – Тучин, 2002;

66 – Тучин, 2010]. Поэтому интервал пребывания на КСО оценивается пятью сутками. К номинальной КСО предъявляется ряд требований. За время пребывания на КСО с борта КА должно быть получено телевизионное изображение предполагаемого района посадки и проверена работоспособность лазерного высотомера автономной системы посадки. В момент начала сеанса посадки, при получении изображения района посадки и проверки работоспособности лазерного высотомера, удаление от поверхности Фобоса не должно превосходить 60 км и должны выполняться условия по освещённости и радиосвязи со станциями слежения в Евпатории и Уссурийске. Если в штатном сеансе работы автономной системы посадка не будет осуществлена и КА останется на прежней орбите, должна быть обеспечена возможность повторения попытки посадки на Фобос.

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

Специалисты Института геохимии и аналитической химии им. М.В.

Вернадского РАН указали несколько интересных для исследования точек поверхности Фобоса. Координаты точек посадки задаются в системе координат, центр которой находится в центре фигуры Фобоса. Положение точек местности определяется сферическими координатами. Широты считаются положительными к северу от экватора. Долгота измеряется в экваториальной плоскости к западу от нулевого меридиана. Сеанс посадки целесообразно начинать с упреждением достижения КА долготы точки посадки. Величина этого упреждения может находиться в диапазоне от 1° до 60°. В момент начала сеанса посадки должны быть обеспечены условия освещённости и радиосвязи с наземными станциями слежения.

По моменту времени прохождения над заданной долготой всегда можно определить момент времени прохождения долготы 270. Поэтому рассмотрим задачу построения такой КСО, которая проходит в заданное время над точкой с долготой 270 и удалена от поверхности на расстояние в диапазоне от 50 до 60 км. Потребуем от этой КСО, чтобы при других прохождениях КА долготы 270 его удаление от начальной точки было минимальным.

Движение по КСО представлено как облёт Фобоса по дрейфующему эллипсу.

Такое описание уже использовано в работах [94 – Henon, 1970; 84 – Benest, 1976;

31 – Коган, 1988; 38 – Лидов, 1994] для получения качественных результатов. В диссертации это представление развивается для получения инструмента для проектирования посадки на Фобос. Большая полуось дрейфующего эллипса ориентирована вдоль орбитального движения Фобоса. КА обегает эллипс в направлении, обратном к орбитальному движению Фобоса, с периодом, близким к периоду обращения Фобоса вокруг Марса. Дрейф эллипса также проходит вдоль оси, ориентированной по орбитальному движению Фобоса.

Для обеспечения успешной посадки нужна такая КСО, дрейф которой был бы минимален. Искомые КСО с минимальным дрейфом обладают тем важным свойством, что прохождение над точкой посадки всегда происходит на одной и той же высоте. Это позволит при необходимости, если процесс посадки не был начат, начать его в одном из следующих прохождений над районом посадки, повторить посадку, произвести телевизионную съёмку района посадки, выполнить измерения высоты до поверхности Фобоса при его облёте.

Вначале рассмотрим численный алгоритм построения КСО, с помощью которого были найдены КСО для удалений на 50 км, 55 км и 60 км. Затем рассмотрим аналитический метод, позволяющий найти начальное приближение.

Этот метод позволил ускорить процесс расчёта множества КСО и вычислить КСО для удалений от поверхности Фобоса в диапазоне от 50 до 200 км.

1.1 Численный алгоритм построения множества КСО с В данном разделе рассмотрен численный алгоритм, позволяющий вычислить три множества КСО, обладающих свойством минимального дрейфа, для удалений от поверхности Фобоса на 50, 55 и 60 км в начальный момент. Каждое такое множество содержит параметры орбиты в зависимости от истинной аномалии Фобоса, которая соответствует моменту времени, в который КА проходит долготу 270°. При поиске КСО в качестве упрощённой модели движения используется плоская эллиптическая задача Хилла в безразмерных переменных. Для ограничения области перебора используются диапазон значений интеграла Якоби и диапазон значений разности угловых скоростей обращения КА вокруг Фобоса и Фобоса вокруг Марса.

1.1.1 Формализация постановки задачи Требуется предложить такой метод выбора начальных условий, который позволил бы для любого заданного момента времени получать КСО, проходящую в этот момент времени на высоте h (50 – 60 км) над точкой поверхности Фобоса с широтой 0° и долготой 270° и обладающую свойством минимального разброса высот относительно h при предыдущих и последующих прохождениях над этой точкой. Задача в такой постановке требует учёта вращения Фобоса вокруг своей оси.

Постановку задачи можно упростить, воспользовавшись тем, что нулевой меридиан Фобоса постоянно обращён к Марсу. Рассмотрим движение КА относительно Фобоса в его орбитальной системе координат. Центр этой системы координат совпадает с центром масс Фобоса. Ось X направлена по линии визирования Марс – Фобос. Ось Y ортогональна оси X, лежит в плоскости орбиты Фобоса и направлена в сторону его орбитального движения. Так как нулевой меридиан Фобоса постоянно обращён в сторону Марса, точкам, расположенным на отрицательной части оси Y, в этой системе координат соответствуют подспутниковые точки с долготой, близкой к 270°. Поэтому будем рассматривать задачу в следующей постановке. Задано положение КА в орбитальной системе координат Фобоса T: X = 0, Y = ( 270 + h ), Z = 0, где 270 — расстояние от центра масс Фобоса до точки экватора Фобоса с долготой 270°. Задавая различные значения вектора скорости, будем получать различные орбиты. При этом только в части случаев КА будет оставаться на КСО [112 – Tuchin, 2003].

Рассмотрим семейство КСО, проходящих через заданную точку T. Для каждой КСО определим множество точек пересечения с отрицательной частью оси y и максимальное удаление точек этого множества от точки T. Будем искать такую КСО, для которой минимально максимальное удаление.

1.1.2 Упрощённая модель движения КА относительно Фобоса Если бы эксцентриситетом орбиты Фобоса (0.015) можно было бы пренебречь, то движение КА относительно Фобоса описывалось бы задачей Хилла [21 – показали, что в этом случае модельные и реальные орбиты сильно различаются используется модель движения, аналогичная задаче Хилла, но с учётом эксцентриситета орбиты Фобоса. В безразмерных переменных уравнения движения имеют вид:

где e – эксцентриситет орбиты Фобоса, – истинная аномалия Фобоса, а точка обозначает производную по.

Безразмерные переменные x и y связаны с положением КА относительно Фобоса X, Y, Z следующими соотношениями:

где pPh – фокальный параметр орбиты Фобоса, Ph – гравитационный параметр Фобоса, M – гравитационный параметр Марса..

Уравнения движения (1.1) могут быть получены из уравнений движения ограниченной задачи трёх тел в форме Нехвила [20 – Дубошин, 1968].

проходящие через точку: x0 = 0, y0 = b0 при заданном 0. В качестве b0 будем рассматривать три значения: 2.456423, 2.654142 и 2.851861, соответствующие удалениям на 50, 55 и 60 км при 0 = / 2. Эти решения будем рассматривать на интервале изменения истинной аномалии от 0 до 0 + 2 104.

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

0.5 < r d, т. к. в рассматриваемом множестве решений r > 0.

Рассмотрим:

Т.к. c0 > d, для 0 применима формула 2.571 п. 4 из [19 – Градштейн, 1971], неполный эллиптический интеграл первого рода:

интеграла.

Используя (1.33), преобразуем (1.32) к виду:

где K ( k ) – полный эллиптический интеграл первого рода.

Таким образом, найдено новое представление этого интеграла через выбранные фазовые переменные a,, A, и свободный параметр.

Используя формулу 8.123 п. 2 [19 – Градштейн, 1971] дифференцирования полного эллиптического интеграла по его модулю и представление производной через полные эллиптические интегралы первого и второго рода, получим:

Далее вычислим производные c0 и d по переменным a, и A :

1.2.7 Усреднённые уравнения движения Система (1.38) имеет интеграл: am Усреднённая система (1.38) совпадает с усреднённой системой (1.29) для нулевого эксцентриситета.

1.2.8 Условие для поиска начального приближения начальных Условие для поиска начального приближения параметров КСО, обладающих искомым свойством, позволяет получить уравнение для m из системы (1.38).

Потребуем, чтобы m =, m = 0. Тогда условие превращается в уравнение, связывающее am и Am :

где Используя (1.35), представим функцию g ( am, Am ) в виде:

Заметим, что при m = и = 2 величину c1 можно было представить полным квадратом: c1 = равным Для определения зависимости a ( A ) нужно решить квадратное уравнение Корни этого уравнения вычисляются по формуле и имеют действительные значения при выполнении условия Расчёты показали, что движение по КСО можно получить только при использовании меньшего корня, т.е. при знаке минус перед радикалом в (1.44).

1.2.9 Алгоритм формирования начального приближения Компоненты вектора положения x0, y0 связаны с параметрами a и A следующим соотношением:

Решая эту систему линейных уравнений при относительно переменных a и A, получим:

При = система (1.45) имеет решение, если выполняется условие:

из решения уравнения: 2a ( A ) A = y0 + 4 x0.

Пусть теперь известно, что КСО должна пройти через точку с координатами x0, y0 в момент времени, в который истинная аномалия Фобоса равна. Тогда, при заданных и, определяем a и A.

Определение начального приближения при поиске начальных условий КСО с искомыми свойствами определяется сканированием по и. Сканирование выполняется в два прохода. На первом проходе сканирование по выполняется на интервале от 0 до 2 с шагом. На втором проходе сканирование выполняется на интервале 0, 0 +, где 0 – точка, в которой была достигнута минимальная ширина кольца на первом проходе. Шаг сканирования на втором проходе составляет на первом и втором проходах одинаковый: от до +, где = 0.3. Шаг сканирования составляет от размера интервала сканирования.

На каждом шаге сканирования определяются a и A. Если A > 0, a ( A ) > 0 и a a ( A ) < a, формируется начальное приближение и вычисляется ширина кольца.

Указанный алгоритм за время расчёта ~15 минут на Pentium 4 позволяет определить начальное приближение с высокой точностью.

1.2.10 Вычислительная процедура поиска минимума Вычислительная процедура поиска минимума построена на основе комбинированного применения метода градиентного спуска и метода золотого сечения [51 – Растригин, 1974]. Сначала применяется градиентный спуск, а затем метод золотого сечения. Алгоритм градиентного спуска дополнен возможностями выхода из локального минимума и идентификации попадания в овраг по одной из переменных.

При использовании градиентного спуска КСО исследуется на небольшом интервале времени, который соответствует интервалу истинной аномалии от 0 до 0 + 2N, где N = 100. Локальный минимум, найденный методом градиентного спуска, уточнялся методом золотого сечения для большего значения N = 10000.

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

Рассмотрим алгоритм градиентного спуска. Как и в [57 – Тучин, 2007] будем рассматривать задачу в следующей постановке. Задано начальное положение КА x0 = 0, y0 = b в момент времени, в который истинная аномалия Фобоса равна 0.

Рассмотрим множество КСО, проходящих через заданную точку, при различных значениях вектора скорости. Для каждой такой КСО определим множество точек пересечения с отрицательной частью оси y и максимальное удаление точек этого множества от начальной точки. Найдём вектор начальной скорости, для которого минимально значение максимального удаления. Численные эксперименты показали, что при использовании градиентного спуска для поиска КСО с минимальным значением целесообразно использовать уравнения в канонических полярных переменных: Q1, Q2, P, P2 :

Начальное положение задается условием:

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

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

где h = 4 105.

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

точка ( P, P20 ) попала в локальный экстремум.

При попадании текущей точки в локальный экстремум вычисляются значения ( P01, P02 + hP 2 ), ( P01, P02 hP 2 ), где hP1 = 4 103, hP 2 = 3 102. Если минимальное значение из этих четырёх значений меньше, чем ( P01, P02 ), происходит переход в соответствующую точку.

Если все четыре значения больше, чем ( P01, P02 ), считается, что текущий локальный экстремум и есть искомый глобальный экстремум, а значения P01, P02 – искомые значения.

( P, P20 ) < ( P h, P20 ), текущая точка попала в овраг по переменной P. В этом случае при переходе к следующей точке должны использоваться только производные по переменной P20, значение P1 обнуляется.

( P, P20 ) < ( P, P20 h ) точка попадает в овраг по переменной P20. В этом случае обнуляется значение P 2. Одновременное попадание в овраг по переменным P и P20 означает, что текущая точка попадает в локальный экстремум. Этот случай уже рассмотрен выше.

Вычисляется величина: m = 2 1 + 2 2. Если эта величина меньше, чем grad = 1012, считается, что достигнут минимум по малому модулю градиента, а значения P01, P02 — искомые.

Дальше ищется значение шага h, при котором Для этого выполняется цикл по величине шага h, начиная со значения, использованного при численном нахождении производных. При каждом новом проходе цикла величина шага h уменьшается вдвое. Если найдено искомое значение шага h, происходит переход к новым значениям обобщённых скоростей P, P20 и повторение действий, начиная с вычисления производных.

Если величина шага h стала меньше, чем 106, считается, что искомый минимум найден, а значения P01, P02 – искомые значения обобщённых скоростей.

Особенность использования метода золотого сечения [51 – Растригин, 1974] состоит в том, что он применяется к функции двух переменных P и P2. Сначала фиксируется P и поиск выполняется по переменной P2. После того, как минимум найден, фиксируется P2 и выполняется поиск по P. Затем выполняется поиск по P и т. д. Другой особенностью применения метода золотого сечения является то, что поиск начинается из точки текущего найденного минимума. Для того чтобы попасть в условия поиска по методу золотого сечения, границы интервала поиска задаются асимметрично относительно текущего минимума. Если C — точка, с которой C r0, C + r0, где r0 — радиус поиска. При поиске минимума методом золотого сечения использованы следующие значения параметров: величины интервалов поиска 0.1 и 0.4 по P1 и P2 соответственно, порог прекращения поиска 0.001 (одинаковый по P1 и P2).

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

Таблица 1.6 – Параметры КСО Начальное удаление Q В таблицах 1.7 – 1.20 приведены значения обобщённых скоростей на начальный момент для начальных удалений от 2.456423 до 8. Значения приводятся с учётом свойства симметрии рассматриваемой системы. Система уравнений (1.10) обладает свойством симметрии: если ( x ( ), y ( ) ) — решение системы (1.10), то ( x ( ), y ( ) ) – тоже решение системы (1.10). В канонических полярных переменных: Q1, Q2, P, P2 симметричные решения уравнения (1.10) Q1( ) ( ), Q2 ) ( ), Q1( При поиске начальных условий рассматривается семейство КСО, проходящее соответствующий значению истинной аномалии 0. Для каждой КСО определяется множество точек пересечения с отрицательной частью оси Y и максимальное удаление точек этого множества от точки T. Для каждого начального значения минимально. Начальные значения обобщённых скоростей обозначим: P01 ( 0 ) и P02 ( 0 ).

Решению уравнения (1.10) с начальными условиями: Q01,, P01, P02 будет Q01,, P01, P02. Очевидно, что максимальные удаления симметричных решений равны между собой. При этом в одном из симметричных решений истинная аномалия растет, в другом убывает. Для периодических решений максимальное удаление от точки T не должно зависеть от направления движения по траектории (вперед или назад). Поэтому, если бы найденные симметричные решения были бы периодическими, максимальные удаления от точки T должны совпадать.

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

P01 ( 0 ) P01 ( 0 ), P02 ( 0 ) P02 ( 0 ), а в таблицах, содержащих результаты расчётов, приведены одинаковые значения P01 и P02 для начальных значений истинной аномалии: 0 и 2 0.

На рис. 1.5 показаны найденные КСО для начальных удалений 3 и 8.

Рисунок 1.5 – Движение по КСО и движение центра эллипса: а) – начальное удаление Q10 = 3, б) – начальное удаление Q10 = Таблица 1.7 – Значения обобщённых скоростей P01 и P02 для начального удаления Q10 = 2. Таблица 1.8 – Значения обобщённых скоростей P01 и P02 для начального удаления Q10 = 2. Таблица 1.9 – Значения обобщённых скоростей P01 и P02 для начального удаления Q10 = 2. Таблица 1.10 – Значения обобщённых скоростей P01 и P02 для начального удаления Q10 = 3. Таблица 1.11 – Значения обобщённых скоростей P01 и P02 для начального удаления Q10 = 3. Таблица 1.12 – Значения обобщённых скоростей P01 и P02 для начального удаления Q10 = 4. Таблица 1.13 – Значения обобщённых скоростей P01 и P02 для начального удаления Q10 = 4. Таблица 1.14 – Значения обобщённых скоростей P01 и P02 для начального удаления Q10 = 5. Таблица 1.15 – Значения обобщённых скоростей P01 и P02 для начального удаления Q10 = 5. Таблица 1.16 – Значения обобщённых скоростей P01 и P02 для начального удаления Q10 = 6. Таблица 1.17 – Значения обобщённых скоростей P01 и P02 для начального удаления Q10 = 6. Таблица 1.18 – Значения обобщённых скоростей P01 и P02 для начального удаления Q10 = 7. Таблица 1.19 – Значения обобщённых скоростей P01 и P02 для начального удаления Q10 = 7. Таблица 1.20 – Значения обобщённых скоростей P01 и P02 для начального удаления Q10 = 8. 1.2.12 Алгоритм проектирования КСО Исходными данными для проектирования являются следующие величины:

h270 – максимальное удаление КА от поверхности Фобоса на долготе 270°, t – момент времени и L – долгота, которую должен проходить КА в этот момент времени. Сначала по максимальному удалению h270 выбирается одна из таблиц 1. – 1.20. Для этого заданное максимальное удаление пересчитывается в безразмерную величину y по формуле:

где – заданное максимальное удаление КА в диапазоне от 50 до 200 км;

h r270 – расстояние от центра Фобоса до точки на его экваторе с долготой 270°, В ходе выполнения работ по баллистико-навигационному обеспечению полёта КА «Фобос-Грунт» по измерениям, полученным на орбите наблюдения, будут уточнены эфемериды Фобоса и его гравитационный параметр. При проектировании КСО будут использованы уточнённые значения Ph и pPh.

По найденной безразмерной величине y и значениям Q10 максимальных удалений выбирается таблица. Выбор выполняется из условия близости y к Q10.

Таблицы построены так, что диапазон удалений от 50 до 60 км, номинального случая в проекте «Фобос-Грунт», покрыт с шагом 5 км, а остальной диапазон от до 200 км, предназначенный для работы при отклонениях от номинального случая, – с шагом 12.5 км.

По заданному моменту времени t определяется t – истинная аномалия Фобоса на этот момент времени. Эта величина используется для выбора строки со значениями обобщённых скоростей из выбранной таблицы. Для этого по начальным условиям, соответствующим строке таблицы, определим долготу ( t ), которую проходит КА при истинной аномалии t. Следует отметить, что каждой строке таблицы соответствует пара начальных условий. Выбирается такая строка таблицы, начальные условия.

По выбранным начальным условиям вычисляется вектор состояния КА в канонических полярных переменных на истинную аномалию t : Q1, Q2, P, P2. На момент времени t по полной модели движения Фобоса в инерциальной системе оскулирующие элементы его орбиты.

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

Проектные расчёты [67 – Тучин, 2002] выполняются в системе координат, центр которой совпадает с центром масс системы Марс – Фобос, плоскость XY совпадает с плоскостью среднего экватора Земли эпохи J2000, ось X направлена в точку весеннего равноденствия Земли этой эпохи, ось Z направлена ортогонально плоскости XY в сторону северного полюса, ось Y дополняет систему координат до правой.

Вектор состояния КА ( rRNB vT ) в орбитальной системе координат Фобоса вычисляется по формулам:

где нулевая матрица 33, единичная матрица 33, переменная, вычисляемая по формуле: t = оскулирующий интеграл площадей орбиты Фобоса, cPh pPh оскулирующий фокальный параметр орбиты Фобоса, оскулирующий эксцентриситет орбиты Фобоса.

ePh Вектор состояния КА в инерциальной системе координат вычисляется по формулам:

где rPh расстояние от центра масс системы Марс – Фобос до центра масс Фобоса на rPh радиальная скорость Фобоса на этот же момент времени.

Матрица CRNB и её производная по времени вычисляются по формулам:

где M PQR – матрица, составленная из векторов Гаусса: P, Q и R.

Векторы Гаусса вычисляются по формулам:

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

Алгоритм проектирования КСО позволяет получить начальные условия, которые обеспечивают прохождение над заданной долготой Фобоса в заданное время. Таблицы начальных условий покрывают интервал от 50 до 200 км.

Глава 2 Баллистика, навигация и управление движением КА на этапе его посадки на поверхность Настоящая глава посвящена проблемам навигации и управления в схеме посадки на Фобос. Результаты, представленные в этой главе, описаны в работах [114 – Tuchin, 2004; 63 – Тучин, 2009].

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

Вернадского РАН. Выполнен анализ схемы посадки и условий ее выполнения.

2.1 Сближение с Фобосом и посадка на его поверхность:

2.1.1 Сближение с Фобосом Для сближения с Фобосом и посадки на его поверхность используются два типа орбит: орбита наблюдения и квазисинхронная орбита (КСО). На первой из орбит выполняются телевизионные наблюдения Фобоса с целью получения точности его эфемерид, которая необходима для выполнения посадки [75 – Шишов, 2008; 76 – Шишов, 2008].

В качестве орбиты наблюдения выбрана орбита со следующими параметрами:

период обращения 8.3 ч; расстояние периария от центра Марса 9905 км; наклонение к плоскости орбиты Фобоса 0°, периодичность сближений с Фобосом – 4 суток;

периодичность сближений с Фобосом, обеспечивающих возможность проведения автономных навигационных измерений – 8 суток; минимальный временной интервал, при котором расстояние от КА до Фобоса менее 1000 км, – 1.7 ч.

После перехода на орбиту наблюдения выполняются следующие операции:

проведение научных экспериментов; выполнение маневров фазирования;

за один месяц до перехода на КСО выполнение измерений по программе, обеспечивающей высокоточный прогноз движения КА относительно Фобоса;

двухимпульсный маневр перехода с орбиты наблюдения на КСО.

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

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

один день подготовки к сеансу посадки;

сеанс посадки или маневр в нештатной ситуации;

повторный сеанс посадки.

2.1.2 Условия посадки Для обеспечения условий посадки КА должен быть приведён в точку, которая находится на высоте не более 60 км над предполагаемым районом посадки. Должны быть выполнены следующие условия при подготовке и на интервале посадки:

получение телевизионного изображения района посадки за несколько суток до начала сеанса посадки;

нахождение угла Солнце – Фобос – КА в диапазоне от 20° до 70° в течение сеанса посадки;

радиовидимость КА со станций слежения в Уссурийске и Медвежьих Озёрах;

остронаправленной антенны (ОНА);

прогноз движения КА относительно Фобоса на момент схода с КСО с ошибками, не превосходящими 3 км по положению и 1 м/с по скорости;

проверка работоспособности основных бортовых систем, обеспечивающих посадку, до её начала;

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

тяговооруженности и запасов топлива КА.

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

бесплатформенный инерциальный блок (БИБ);

звёздный прибор (БОКЗ-МФ);

лазерный высотомер (ЛВ);

доплеровский измеритель скорости и дальности (ДИСД);

телевизионная система.

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

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

Лазерный высотомер обеспечивает измерения расстояния до поверхности по четырем лучам. Он выполняет навигационные измерения на участке от момента схода с КСО до высоты 500 м. Доплеровский измеритель скорости и дальности выполняет измерения дальности до подстилающей поверхности в направлении четырёх лучей и проекции вектора скорости КА относительно Фобоса на направления этих лучей. Доплеровская система работает с высоты 3 км. Измерения лазерного высотомера и доплеровской системы используются навигационной задачей, входящей в бортовой комплекс управления посадкой, для определения вектора состояния КА, которым корректируется вектор состояния, получаемый методом счисления пути. Измерения дальностей по лучам, получаемые ЛВ и ДИСД, используются также для определения нормали к подстилающей поверхности.

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

2.1.4 Схема управляемой посадки Схема посадки на Фобос (рис. 2.1) потребовала анализа траектории с использованием четырёх участков:

схода с КСО;

перелёта с КСО в точку, расположенную над районом посадки;

вертикального спуска;

прецизионного торможения.

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

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

При этом в горизонтальной плоскости накапливается ошибка по положению.

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

Рисунок 2.1 – Схема управляемой посадки в СК, связанной с фигурой Фобоса На участке прецизионного торможения постепенно гасится вертикальная составляющая скорости до величины, с которой допускается соприкосновение с поверхностью: 1.5 – 2 м/c. При этом боковая составляющая вектора скорости не должна превосходить 1 м/с.

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

До начала сеанса посадки выполняется разворот, обеспечивающий заданную ориентацию. Она определяется из условий:

«захвата» поверхности Фобоса лучами ЛВ;

освещённости солнечных батарей в ходе и после посадки;

работы привода ОНА, обеспечивающей связь с Землей.

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

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

по длительности всей операции посадки;

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

по энергетическим затратам.

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

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

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

Моделирование процесса посадки показало, что длительность перелёта от момента схода с КСО до начала участка вертикального спуска составляет ~30 мин.

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

Выбор точки посадки, варианта схода с КСО, высоты приведения над точкой посадки, расчёт импульса схода с КСО выполняется на Земле. Перед выполнением сеанса посадки на борт передаются следующие данные:

момент схода с КСО и длительность перелёта от момента схода до прицельной импульс схода и положение прицельной точки;

заданная ориентация КА, векторы состояния КА и Фобоса на момент схода с моменты возможных коррекций на интервале перелёта.

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

2.1.5 Управление включением двигателей Управление движением КА на участке посадки обеспечивают 20 двигателей.

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

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

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

Заданный импульс и набранная характеристическая скорость преобразуются к «повёрнутой» системе координат. Наличие невязки между компонентой вектора характеристической скорости и соответствующей ей компонентой заданного импульса означает, что необходимо создать ускорение вдоль соответствующей оси «повёрнутой» системе координат. Это ускорение создаётся включением определённых двигателей.

Заданная и текущая ориентации также рассматриваются в «повёрнутой»

системе координат. Ориентация определяется углами Эйлера. Наличие невязки между текущим и заданным значениями углов означает, что КА нужно повернуть относительно соответствующей оси.

Таким образом, управление движением КА по каждому каналу реализует обратную связь по невязке. Закон управления представляет собой релейное управление с гистерезисом (чтобы избежать скользящего режима). Такой закон определяется двумя порогами по каждому направлению. Больший по абсолютной величине порог ср – порог на включение (срабатывание), меньший откл – на отключение (рис. 2.2).

Рисунок 2.2 – Вид нелинейности в цепи обратной связи 2.2 Алгоритмы управления движением КА 2.2.1 Бортовые алгоритмы навигации и управления Основными бортовыми алгоритмами навигации и управления являются алгоритмы счисления пути, определяющие вектор состояния КА, его ориентацию и вектор угловой скорости, алгоритмы навигационной задачи и алгоритм расчёта коррекции.

Функциональная схема работы алгоритмов навигационного обеспечения управления посадкой показана на рис. 2.3. Вектор состояния КА и соответствующая ему ковариационная матрица на начало сеанса посадки рассчитываются на Земле по данным наземных и бортовых траекторных измерений. Эти данные передаются на борт и являются исходными для алгоритма счисления пути.

Алгоритм счисления пути по измерениям векторов кажущихся ускорений интегрированием уравнений движения КА в ограниченной задаче трёх тел (Марс, Фобос, КА) определяет оценку вектора состояния КА на текущий момент времени.

Эта оценка заведомо будет включать ошибку по положению и скорости КА относительно Фобоса, содержащуюся в начальных условиях движения КА, полученных с Земли.

Исходный вектор состояния и ковариационная матрица Алгоритм управления Рисунок 2.3 – Навигационное обеспечение управления посадкой По мере сближения КА с Фобосом от лазерного высотомера и доплеровской системы поступает более точная информация о положении и скорости КА относительно Фобоса. По этой информации решается навигационная задача и корректируется оценка вектора состояния КА, полученная методом счисления пути:

алгоритм счисления пути получает новые начальные условия.

2.2.2 Уравнения движения КА относительно Фобоса Алгоритм решения навигационной задачи и алгоритм расчёта коррекций используют уравнения движения КА относительно Фобоса. Уравнения записаны в его орбитальной СК, центр которой совпадает с центром масс Фобоса. Ось Ox направлена по линии визирования Марс – Фобос. Ось Ox2 ортогональна оси Ox1, лежит в плоскости орбиты Фобоса и направлена в сторону его орбитального движения. Ось Ox3 дополняет систему до правой. Орбитальная СК Фобоса ещё удобна и тем, что переход из неё в СК, связанную с фигурой Фобоса, выполняется линейным преобразованием, коэффициенты которого постоянны на интервале посадки. Коэффициенты этого линейного преобразования вычисляются на Земле при подготовке сеанса посадки Далее будем использовать обозначения, введённые в первой главе. Для орбиты Фобоса: e – эксцентриситет, pPh – фокальный параметр, cPh – интеграл площадей, – истинная аномалия. Гравитационные параметры Марса и Фобоса обозначим M, Ph. Для КА введём два вектора: x = ( x1, x2, x3 ) – вектор положения a = ( a1, a2, a3 ) – вектор кажущегося ускорения в орбитальной СК Фобоса. Используя относительно Фобоса в его орбитальной СК и нормированных переменных:

Векторы положения x RF и скорости x RF КА в орбитальной СК связаны с векторами положения и скорости в нормированных переменных следующими соотношениями:

Методические ошибки модели движения (2.1), (2.2) на участке схода длительностью ~40 минут не превышают 100 м по положению и 0.05 м/с по скорости.

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

где матрица A() вычисляется вдоль решения системы (2.1), E 3 и 03 – единичная и нулевая квадратные матрицы размерности три.

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

Матрица частных производных текущего вектора состояния КА в орбитальной СК Фобоса по начальному вектору состояния в этой же системе координат вычисляется по формуле:

где D ( t ) – матрица частных производных вектора состояния в орбитальной СК на момент времени по компонентам вектора состояния в нормированных переменных. В формуле (2.4) моменту времени t соответствует значение истинной аномалии, а моменту времени t0 – значение 0. Матрицы D ( t ) и D 1 ( t0 ) вычисляются по формулам:

2.2.3 Расчёт коррекций на участке перелёта от момента схода с КСО до точки начала вертикального спуска Матрица (2.4) позволяет выполнить расчёт коррекции на участке перелёта после схода с КСО в точку начала участка вертикального спуска. Если моменту времени t0 соответствует момент проведения коррекции, а моменту времени t момент выхода в точку начала вертикального спуска, то матрица (2.4) позволяет связать в линейном приближении невязки по компонентам положения на момент t с поправками к вектору скорости на момент времени t0, которыми компенсируются невязки по положению. Расчёт импульса коррекции производится итерационно. На каждом шаге итерационного процесса уточняется вектор приращения скорости, обеспечивающий приведение КА в точку, заданную в орбитальной СК Фобоса. Для расчёта коррекции необходимо не более трёх итераций. Первая итерация позволяет определить вектор приращения скорости с точностью до единиц сантиметров в секунду, а третья – с точностью до долей миллиметров в секунду.

2.2.4 Определение вектора состояния КА по измерениям лазерного высотомера и доплеровской системы Навигационная задача (определение вектора состояния КА по измерениям лазерного высотомера и доплеровской системы) решается классическим методом наименьших квадратов [12 – Аким, 1963]. Измеряемыми функциями являются расстояние до поверхности Фобоса вдоль заданного направления и проекция вектора скорости на заданное направление.

Введём обозначения: x – искомый вектор состояния КА на момент времени t0 ; x 0 – априорный вектор состояния КА, полученный методом счисления пути на момент времени t0 ; C0 – ковариационная матрица ошибок оценки вектора состояния x 0 ; N – число измерений; wi вес измерения; imeas измеренное значение расстояния до поверхности в заданном направлении imeas или проекции скорости КА относительно Фобоса на заданное направление imeas ; icalc расчётный аналог измерения: icalc или icalc. Алгоритм решения навигационной задачи основан на минимизации функции, представляющей собой сумму квадратов взвешенных невязок расчётных и измеренных значений:

Поиск минимума функции приводит к решению нелинейной системы уравнений методом Ньютона, который сводится к серии последовательных прибавляемых к оценке вектора состояния x ( приближении, решается система нормальных уравнений:

Матрица и правые части этой системы формируются для приближения с номером s по следующим формулам:

где x ( измеренного значения по вектору состояния КА в СК, связанной с фигурой Фобоса;

вектору положения на момент t0, вычисляемая по формуле (2.4). Матрица M RF_Ph ( ti ) для пересчета вектора положения из орбитальной СК Фобоса в СК, связанную с его фигурой, имеет простой вид: M RF_Ph ( ti ) = diag(1, 1, 1) + G, где G матрица поправок, которая задаётся постоянной на участке посадки.

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

Расчётный аналог измеренного значения проекции вектора скорости на заданное направление вычисляется по формуле:

где k i – единичный вектор заданного направления луча в орбитальной СК Фобоса;

ti – момент времени измерения i.

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

Для поверхности Фобоса используются две аппроксимации: референцэллипсоид (РЭ) с длинами полуосей a, b, c и цифровая поверхность Фобоса (ЦПФ), которая задаётся таблицей с шагом по широте и долготе. Для каждого узла дана величина радиус-вектора от центра Фобоса до его поверхности.

Основная идея алгоритма построена на представлении ЦПФ в виде отклонений от поверхности РЭ в направлении нормали к ней.

Пусть re = ( xe, ye, ze ) – точка, лежащая на поверхности РЭ. Тогда уравнение нормали к поверхности: r = re + pn, n = e, 2, e, p – коэффициент. Пусть теперь rs = ( xs, ys, zs ) – точка на ЦПФ. Точки re и rs лежат на нормали к поверхности РЭ.

Требуется получить координаты точки re на поверхности РЭ. Из параметрического выражения в уравнение эллипсоида, получим уравнение:

Это уравнение нужно разрешить относительно параметра p.

В алгоритме, который будет изложен ниже, несколько раз производится переход от точки на РЭ к точке на ЦПФ, причем обе точки лежат на нормали к РЭ.

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

1. После нахождения матрицы P вычисляются узловые точки на поверхности РЭ.

Их угловые координаты e и e распределены неравномерно (в отличие от угловых координат ЦПФ). Так как в реальном времени для вычисления промежуточных значений коэффициента p будет использоваться линейная интерполяция, нужно вычислить новую матрицу Pun коэффициентов для перехода от точек РЭ к точкам ЦПФ в равномерной сетке угловых координат.

Для этой цели использована так называемая триангуляция Делоне [55 – Сначала изложим основную идею алгоритма. Пересечение луча с ЦПФ ищется итерационно. Первая точка rb,1 на луче определяется как пересечение луча с РЭ.

Вторая и последующие точки rb,i +1, i 1 вычисляются при выполнении следующих действий:

на ЦПФ определяется точка rs,i, лежащая на нормали к РЭ, проходящей через точку rb,i ;

через найденную точку rs,i проводится плоскость, ортогональная нормали к РЭ;

находится точка rb,i+1 пересечения луча с этой плоскостью.

Приведём вспомогательный алгоритм нахождения точки на ЦПФ, лежащей на нормали к точке на РЭ. Пусть re = ( xe, ye, ze ) – координаты точки, лежащей на РЭ.

Вычисляем для re угловые координаты ( e, e ). По матрице Pun вычисляем для узла ( e, e ) коэффициент перехода p. Для этого используем линейную интерполяцию.

Тогда rs = re + pn.

Теперь мы готовы изложить итеративный алгоритм нахождения точки пересечения луча ЛВ с ЦПФ (рис. 2.3). Исходные данные: r0 – положение КА в системе координат, связанной с Фобосом, k – направляющий вектор лазерного луча в этой же системе координат. Требуется найти rs – точку пересечения поверхности Фобоса и лазерного луча.

Рисунок 2.3 – Нахождение точки пересечения лазерного луча и поверхности Фобоса 1-й шаг. Вычисление rb,1 – точки пересечения луча и РЭ. Точка rb,1 является первым приближением к точке rs. Здесь может возникнуть особая ситуация:

пересечение с эллипсоидом не найдено. Тогда алгоритм аварийно завершает работу.

2-й шаг. Обозначим rs,1 – точку ЦПФ, лежащую на нормали к РЭ, проходящей через точку rb,1. По вспомогательному алгоритму найдём вектор d, соединяющий rb,1 с точкой rs,1. Через rs,1 проведём плоскость, ортогональную вектору нормали. Эта плоскость пересекает луч в точке rb,2 = rb,1 + k, которая является вторым приближением к искомой точке.

3-й и последующие шаги. Зная точку rb,i, i 2, лежащую на лазерном луче, найдём её ортогональную проекцию на РЭ – точку re,i. Для этого решаем уравнение (5) относительно p итерационным численным методом. По вспомогательному алгоритму отыскиваем точку rs,i, лежащую на ЦПФ. Проведём через rs,i плоскость, ортогональную нормали к эллипсоиду в точке re,i. Пересечение этой плоскости с лазерным лучом даёт точку rb,i+1 = r0 + k – очередное i + 1-е приближение к искомой точке.

Численные эксперименты показали, что на точке rb,4 итерации можно остановить, т.к. последующие точки лишь незначительно от неё отличаются.

Частные производные дальности от r0 до rb,i по компонентам вектора r приближённо находятся как частные производные дальности от r0 до re,1 по компонентам вектора r0, которые, в свою очередь, вычисляются по аналитическим формулам.

Матрица Pun представляется в памяти бортового компьютера в компактном виде. Так как разброс её элементов невелик (от –13.6 до 18.2), они могут быть представлены в формате short int. Общий объём памяти, занимаемой Pun, составляет 32942 байта.

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

Лавочкина при непосредственном участии автора диссертации в системе управления КА «Фобос-Грунт», разрабатываемой в НПО им. С.А. Лавочкина. Эти алгоритмы описаны в научно-технических отчётах [131 – НТО 5-07-05; 133 – НТО 5-036-07;

134 – НТО 5-11-04; 136 – НТО 5-12-06; 137 – НТО 5-026-07; 138 – НТО 5-016-08;

140 – НТО 5-11-05].

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

Основная задача работает в рамках такта управления КА. Фоновая задача выполняет расчёты, распределённые по тактам управления.

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

Фоновая задача определяет вектор состояния КА по измерениям лазерного высотомера и доплеровской системы, выполняет расчёт коррекций на участке перелёта от момента схода до начала участка вертикального спуска.

Тестирование программного обеспечения управления посадкой выполнялось на стенде ИПМ им М.В. Келдыша. Стенд содержит макет бортовой машины и персональный компьютер для моделирования внешней среды и систем КА. В систему управления на макете бортовой машины поступают измерения, получаемые моделями измерительных приборов: БИБ, ЛВ, ДИСД, БОКЗ-МФ. Система управления определяет номера включаемых двигателей. В персональном компьютере моделируется работа двигательной установки с учётом переходных процессов и определяется векторы тяги и момента, создаваемые двигателями. Далее интегрируются уравнения движения центра масс КА, Фобоса и уравнения движения КА вокруг центра масс. Полученные параметры используются для моделирования выхода измерительных приборов.

Схема работы системы моделирования сближения и посадки КА на Фобос показана на рис. 2.5. Система состоит из блоков, моделирующих работу систем КА, и блоков, моделирующих внешнюю среду.

В ходе моделирования процесса посадки КА на Фобос моделируется работа следующих систем КА: измерительных систем (БИБ, ЛВ, ДИСД, БОКЗ-МФ), привода остронаправленной антенны (ОНА), двигательной установки (ДУ).

Модель внешней среды содержит блоки:

интегрирования движения центра масс (ЦМ) КА и вокруг ЦМ;

модели движения ЦМ Фобоса и его движение вокруг ЦМ;

модели рельефа Фобоса.

Система управления посадкой, реализуемая на макете бортовой машины, получает информацию от моделей БИБ, ЛВ и ДИСД. От модели БИБ поступают приращения углов поворота относительно осей приборной системы координат (ПСК) и приращения кажущейся скорости вдоль осей ПСК за интервал времени между опросами. Система управления посадкой восстанавливает ориентацию осей связанной СК относительно J2000 с использованием начальной ориентации, определённой звёздным прибором БОКЗ-МФ. От модели ЛВ поступают измерения дальности по лучам ЛВ в приборной СК. От модели ДИСД в систему управления поступают измерения дальности и скорости по лучам ДИСД в приборной СК.

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

Модель ДУ получает на вход шкалу из 20 бит, каждый бит которой соответствует двигателю: 1 означает, что двигатель должен быть включен, а 0 – выключен. В результате работы модели ДУ определяются ускорения и моменты в связанной СК.

Блок интегрирования уравнений движения ЦМ КА и его движения вокруг центра масс получает ускорения и моменты, создаваемые ДУ, как входную информацию. Выходом этого блока являются вектор состояния и вектор кажущихся ускорений в марсоцентрической СК J2000, ориентация связанных осей КА и угловые скорости. Эта информация поступает на вход моделей измерительных систем: БИБ, ЛВ и ДИСД.

Блок, моделирующий рельеф Фобоса, обеспечивает расчёт расстояния до центра фигуры Фобоса как функцию широты и долготы. Эта информация используется моделями ЛВ и ДИСД.

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

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

Рисунок 2.5 – Схема работы моделирования сближения с Фобосом и посадки на его поверхность Глава 3 Определение параметров движения КА по результатам измерений при наличии немоделируемых ускорений Задача определения параметров движения КА является одной из основных задач, решаемых в ходе управления его полётом. При решении этой задачи часто возникает ситуация, в которой определение параметров движения КА надо выполнять на фоне работы двигателей. В качестве примера можно привести следующие задачи:

контроль участка выведения КА на орбиту искусственного спутника Земли [79 – Akim, 1999; 10 – Аким, 1999; 28 – Тучин, 1991; 69 – Тучин, 1993; 70 – Тучин, 1991];

оперативная оценка исполнения импульсов и прогноза падения орбитального комплекса по измерениям наземных средств на фоне работы двигательной установки [125 – НТО 5-01-02];

определение параметров движения КА с электроракетной двигательной установкой (ЭРДУ) [68 – Тучин, 2010].

К вопросам, которые рассматриваются в диссертации, в первую очередь, относятся проблемы решения задач навигации и управления полётом с включённой ЭРДУ. Фактическое ускорение, создаваемое ЭРДУ, отличается от модели этого ускорения, заложенного в расчёты. Имеются ошибки величины и ориентации вектора тяги ЭРДУ в пространстве.

динамических систем, в которых помехи имеются не только в измерениях, но и влияют на поведение самого объекта. Такие модели исследованы в основном в рамках линейных моделей в общей теории систем [13 – Аоки, 1967; 35 – Красовский, 1968; 43 – Льюнг, 1991; 47 – Острём, 1973; 36 – Ли, 1966; 27 – Калман, 1971; 93 – Haupt, 1995] Применение этих моделей и методов в задачах определения движения КА требует использования соответствующих нелинейных моделей и учёта особенностей уравнений динамики и измеряемых функций.

Результаты, представленные в данной главе, описаны в работе [62 – Тучин, 2009] в части алгоритмов определения параметров движения КА и в работе [58 – Тучин, 2004].

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

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

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

метод наименьших квадратов;

расширенный фильтр Калмана.

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

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

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

Кроме того, метод требует, чтобы не было длительных интервалов времени, в которых нет измерений.

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

Минимизация функционала выполняется итерационно. На каждом шаге итерационного процесса определяется поправка к искомым параметрам.

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

Во втором разделе приведён алгоритм оценки вектора состояния в случае отсутствия возмущений.

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

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

Если в качестве оцениваемого вектора состояния выбрать вектор состояния на конец мерной базы, то неучтённый шум будет приводить к увеличивающимся ошибкам модели по мере перемещения от конца мерной базы к её началу. Суть метода мешающих параметров состоит в учёте этой нарастающей ошибки модели в весовой матрице измерений. Алгоритм оценки вектора состояния с использованием метода мешающих параметров рассмотрен в четвёртом разделе главы. Рассмотрены варианты алгоритма для двух типов возмущений: белого шума и случайных величин, постоянных на всем интервале. Алгоритм оценки вектора состояния с использованием метода мешающих параметров целесообразно применять при решении задач оценки точности определения параметров движения КА. Этот алгоритм позволяет оценить воздействие шума при приближённых представлениях о его статистических характеристиках.

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

3.1 Алгоритм оценки вектора состояния и суммарных воздействий возмущений между измерениями 3.1.1 Постановка задачи Представим модель движения КА на интервале t0, t N в виде суммы опорного движения и движения относительно опорного:

где x D – вектор состояния опорного детерминированного движения;

x p – вектор состояния движения относительно опорного.

Опорное движение и движение относительно опорного описываются следующими системами дифференциальных уравнений:

где F ( t, x ) – вектор-функция;

B(t ) – матрица, описывающая воздействие шума на систему;

w (t ) – случайный процесс с независимыми приращениями.

Пусть белый шум (t ) с матрицей интенсивности Q ( t ) является производной случайного процесса w ( t ).

Начальные условия для (3.1) задаются априорным вектором x 0 и его ковариационной матрицей P0.

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

каждого момента времени ti справедливо где i – случайный вектор, имеющий нулевое математическое ожидание и ковариационную матрицу Ri.

Запись в качестве параметра x ( ) функции i означает, что функция i зависит не от мгновенного значения вектора состояния, а от функции x ( t ), которая может быть представлена в виде (3.1), (3.2).

Рассмотрим задачу в линейной постановке.

3.1.2 Линейный случай В линейной постановке модель движения (3.1), (3.2) примет вид линейного стохастического дифференциального уравнения где A(t ) – квадратная матрица порядка n n, элементы которой являются непрерывными функциями времени t.

линейными функциями вектора состояния x ( t ). В каждый момент времени ti справедливо соотношение где – вектор параметров размерности ri, измеряемых в момент времени ti ;

H i ( ti ) – матрица размерности ri n ;

– последовательность независимых случайных векторов, имеющих нулевое математическое ожидание и ковариационную матрицу Ri .

Зависимость между векторами состояния системы в дискретные моменты времени ti, определяемая дифференциальным уравнением (3.4), может быть выражена при помощи его разностного аналога, определяемого соотношениями Здесь – фундаментальная матрица, удовлетворяющая матричному уравнению при начальном условии ( ti, ti ) = E, где E — единичная матрица размерности n.

{v ( ti ), i = 0,1,2,..., N } – последовательность случайных векторов с нулевым математическим ожиданием и ковариационными матрицами Qi, вычисляемыми по формуле Случайный вектор v ( ti ), i = 0,1,2,..., N связан с шумом (t ) соотношением Этот случайный вектор v ( ti ) интерпретируется, как суммарное воздействие возмущений на интервале индексами i и i + 1.

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

Требуется получить оценку вектора состояния x i дискретной динамической системы, которая описывается следующим соотношением при априорно заданной оценке начального вектора состояния x 0 и ковариационной матрице этой оценки P0. Измеряемые векторы zi, i = 1,..., N,... связаны с векторами состояния уравнениями где i — случайный вектор размерности ri с нулевым математическим ожиданием и ковариационной матрицей R i.

Построим оценку вектора состояния и суммарных возмущающих воздействий дискретной динамической системы (3.11, 3.12) по мерной базе, содержащей N измерений. Используем метод наименьших квадратов. Критерием качества оценки является квадратичная форма следующего вида где x i, N – оценка вектора состояния на момент ti, i = 0,..., N с использованием N vi, N – оценка вектора суммарных возмущений vi при использовании N Квадратичная форма содержит члены трёх типов:

квадраты невязок измеренных и расчётных значений, отнесённые к априорно известным среднеквадратическим отклонениям ошибок измерений;

отнесённые к априорно известным средним значениям суммарных возмущений на тех же интервалах;

квадрат взвешенного отклонения априорно заданного вектора состояния от его расчётного значения.