Известия РАН. Механика твердого тела, 2023, № 1, стр. 68-75

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

Один из основных вопросов при решении навигационных задач – какой математический аппарат использовать для параметризации углового движения [8]? Как правило, это матрицы или кватернионы. Оба подхода имеют свои преимущества и недостатки. В процессе навигации кватернионы или матрицы используются в двух основных операциях: перепроектирование, то есть перевод вектора из одной системы координат в другую, и интегрирование кинематического уравнения, которое показывает как матрица или кватернион перехода изменяется со временем. Матрицу перехода можно параметризовать самолетными углами: крен, тангаж, рыскание. При этом кинематические уравнения для углов нелинейны и могут вырождаться. Непосредственное интегрирование кинематического уравнения для матрицы перехода не содержит эту проблему, однако возникает вопрос об отслеживании ортонормированности матрицы. Кватернионное кинематическое уравнение является линейным и не вырождается, но при этом необходимо следить за нормой. В вопросах выставки и коррекции БИНС предпочтительно использовать те параметры, на которые не накладываются условия нормировки, а также уравнение перепроектирования которых линейно. Несмотря на то, что перепроектирование с помощью кватернионов этим условиям не удовлетворяет, возникающих осложнений можно избежать. В данной работе показано как это можно сделать с помощью алгебраических свойств кватернионов. Как альтернатива, существует описание задачи об определении ориентации аппарата на спуске с использованием матрицы перехода [9]. Однако, оставлен без внимания вопрос: как избежать осложнений, возникающих из-за условий нормировки, налагаемых на компоненты матрицы?

Общая навигационная задача ставится следующим образом: определить все составляющие вектора состояния (ВС) – координаты, скорость, кватернион ориентации – в данный момент времени. Традиционно эта задача решается численным решением системы дифференциальных уравнений, а также использованием специальных систем уравнений – фильтров [10]. Они возникли на стыке теории управления и теории вероятностей и впервые описаны в статьях Р. Калмана в 1960 годах. Термин “фильтр” обозначает систему дифференциальных или разностных уравнений для ВС, в которую входит информация от внешних измерителей о какой-либо части ВС. В данном случае ВС содержит вектор координат, вектор скорости и кватернион ориентации. Внешней информацией служит скорость АСН. Использование фильтров в случае системы нелинейных дифференциальных уравнений затруднено, поскольку изначально они предназначались для линейных систем. Для нелинейных систем возникает необходимость применения дополнительных итерационных методов, позволяющих искать минимум квадратичного критерия, например метод градиентного спуска. Более того, фильтрация подразумевает точное знание параметров шума измерителей и наблюдаемой системы, что трудно осуществить на практике. Поэтому в данной работе рассматривается наиболее классический метод наименьших квадратов (МНК), применение которого не требует знания свойств погрешности приборов. При этом воздействие систематических составляющих должно быть минимальным, либо учитываться отдельно. Одна из причин, по которой был выбран именно МНК, заключается в следующем принципе решения практических вопросов управления/навигации: необходимо свести задачу к наиболее изученной математической форме, имеющей единственное решение. Такой математической формой является система линейных алгебраических уравнений (СЛАУ). Практика показывает, что если МНК применим, то он сводится именно к СЛАУ, что продемонстрировано в данной работе.

2. Постановка задачи. Система дифференциальных уравнений (СДУ) БИНС выглядит следующим образом

Здесь ${\mathbf{a}},{\mathbf{\omega }}$ – негравитационное ускорение и угловая скорость, измеренные акселерометром и ВОГ. Основная часть ВС ${\mathbf{r}},{\mathbf{v}},R$ – координаты скорость в замороженной гринвичской системе координат (ЗГСК) и кватернион перехода из ЗГСК в связанную с аппаратом систему координат (ССК). Вспомогательная часть ВС ${\mathbf{y}}$ – интеграл негравитационного ускорения в ЗГСК. Используемые обозначения: кватернионное умножение $ \circ $, сопряженный кватернион $\tilde {R}$, скалярное и векторное произведение (,), [,]. Понятие “негравитационное ускорение” используется вместо традиционного “кажущееся ускорение” для того, чтобы подчеркнуть суть измерений акселерометра.

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

с начальным условием ${{R}_{0}}$, отличается от решения того же уравнения с начальным условием $K \circ {{R}_{0}}$ умножением на кватернион $K$ слева. Можно убедиться непосредственной подстановкой

В связи с этим корректирующий кватернион K считается постоянным и в вектор состояния не входит в данной постановке. Полная информация об ориентации аппарата содержится в произведении $K \circ R$ после того, как найден K. Более того, точность начальной ориентации ${{R}_{0}}$ не влияет на конечный результат, если удается найти кватернион K. Вектор гравитации g(r) моделируется с учетом второй зональной гармоники разложения гравитационного потенциала по сферическим функциям. В качестве инерциальной системы координат, в которой интегрируется СДУ БИНС, используется ЗГСК.

Решение системы дифференциальных уравнений БИНС численным методом, например Рунге-Кутты четвертого порядка, позволяет осуществить основную навигационную задачу – определение ВС в данный момент времени. Задание начальных условий интегрирования, то есть ${{{\mathbf{r}}}_{0}},{{{\mathbf{v}}}_{0}},{{R}_{0}}$ в начальный момент времени – это задача выставки БИНС. Поскольку численное интегрирование, а также показания приборов не являются идеально точными, со временем ВС становится известен с некоторой погрешностью $\Delta {\mathbf{r}},\Delta {\mathbf{v}},\Delta R$. Задача коррекции заключается в устранении этой погрешности. Показания АСН с высокой точностью определяют непосредственно координаты и скорость, что по сути является коррекцией и выставкой r, v. Кватернион ориентации R по показаниям АСН определяется с низкой точностью и только при определенном наведении антенн на навигационные спутники, а это невозможно гарантировать во время спуска. Например, в [9] указано, что точность ориентации после накопления измерений в течение 5 минут составляет 2 градуса, согласно моделированию и летным испытаниям аппаратуры АСН-К. Поэтому в СДУ БИНС вводится корректирующий кватернион K, который необходимо определить косвенно.

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

Для интегрирования гравитационного ускорения по времени используется линейная интерполяция для концов интервала $[{{t}_{i}},{{t}_{{i + 1}}}]$, поскольку только в эти моменты с высокой точностью известны ${{{\mathbf{r}}}_{i}} = {\mathbf{r}}\left( {{{t}_{i}}} \right),{{{\mathbf{r}}}_{{i + 1}}} = {\mathbf{r}}\left( {{{t}_{{i + 1}}}} \right)$:

Данный подход к вычислению интеграла гравитационного ускорения ограничивает применимость в том смысле, что точность ухудшается по мере увеличения интервала $[{{t}_{i}},{{t}_{{i + 1}}}]$. Интеграл скорости:

Интеграл негравитационного ускорения можно рассчитать по соответствующей переменной ВС:

Таким образом, если совершить замену:

можно прийти к формуле, в которой неизвестен только кватернион K:

В [8, 11] указан способ определения кватерниона перехода по двум неколлинеарным векторам, заданным в своих базисах. Основная проблема данного подхода заключается в том, что на протяжении небольшого промежутка времени векторы ${{{\mathbf{u}}}_{i}}$ и ${{{\mathbf{p}}}_{i}}$ не меняют направления. Это наталкивает на использование не двух векторов, а всех накопившихся в полете. Тогда формально составленный функционал для $N$ измерений выглядит следующим образом:

Здесь новая независимая переменная ${{\lambda }}$ – множитель Лагранжа. Он необходим, поскольку среди четырех компонент кватерниона $K = {{k}_{0}} + {\mathbf{k}}$ независимых только три в силу условия нормировки:

То есть возникает задача поиска условного экстремума. Ее решение в режиме реального времени не представляется возможным, поскольку она сводится к нелинейной системе из 5 алгебраических уравнений. Для того чтобы избавиться от уравнения связи и от нелинейности выражение

умножается справа на $K$, а также делится на ${{k}_{0}}$. В результате:

Поскольку требуется, чтобы все

необходимо минимизировать сумму соотношений:

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

Здесь a, b – некоторые постоянные векторы. Итак, произведение сопряженных кватернионов:

И производные отдельных компонент:

После упрощений:

Здесь E – единичная матрица. Итоговое условие стационарности точки e у функции $F({\mathbf{e}})$ сводится к СЛAУ:

Определив из этой системы вектор ${\mathbf{e}}$, можно вычислить кватернион $K$ по формуле:

В случае, если СЛАУ не имеет решения, то есть ее определитель равен 0, можно принять ${{k}_{0}} = 1$ и ${\mathbf{k}} = 0$, то есть коррекция отсутствует. Для того, чтобы стационарная точка была минимумом необходимо, чтобы вторая производная:

была положительно определенной матрицей. Для доказательства можно проверить свойство соответствующей квадратичной формы: для любых ${{{\mathbf{u}}}_{i}},{{{\mathbf{p}}}_{i}}$. Аналогично для всей суммы. В связи с этим, единственная стационарная точка F – это минимум.

В доказательство того, что приведенный алгоритм согласуется с уже существующей теорией, можно проверить соответствие с формулой определения кватерниона ориентации по двум векторам [8, 11]:

для случая N = 2. Рассмотрение производных отдельных членов суммы до раскрытия: с подстановкой ${\mathbf{e}}$, убеждает в том, что формулы взаимно однозначны. Для i = 1, 2: в силу того, что:

Таким образом

4. Анализ результатов. Приводимый алгоритм определения углового положения КА имеет ряд преимуществ и недостатков по сравнению с аналогичными методами. Основное преимущество – независимость от свойств шумов приборов и входных воздействий. В противоположность фильтру Калмана, в котором ковариационные матрицы используются для интегрирования дифференциального уравнения матрицы ошибок. Более того, размер матрицы ошибок в задаче с ВС ${\mathbf{r}},{\mathbf{v}},R$, состоящим из 10 компонент – 10х10, ее интегрирование в реальном времени может потребовать значительных вычислительных затрат. Еще одно преимущество данного метода – четкий критерий корректности выполнения, а именно определитель СЛАУ не равен 0. Также использование МНК в данном случае позволяет хранить сумму измерений, а не каждое измерение в отдельности, и модифицировать ее при поступлении новых данных.

Недостаток метода коротко можно обозначить так: МНК находит среднее значение K. Это означает, что если ошибка определения кватерниона ориентации $\Delta R(t)$ меняется со временем, то МНК позволит устранить не текущую ошибку, а именно среднюю $\overline {\Delta R(t)} $ за некоторый временной интервал сбора данных T. Например, при длительной работе алгоритма $\Delta R(t)$ медленно меняется в связи с ошибками численного метода решения СДУ БИНС, а также из-за погрешности измерений ${\mathbf{\omega }}$. Также может возникнуть сильное кратковременное изменение $\Delta R(t)$ из-за выхода ВОГ за линейный диапазон измерений и в этом случает ошибка также не может быть полностью устранена. Поэтому МНК не является легко применимым автономным методом, для его длительного использования требуются дополнительные настройки типа обнуления данных или весовых коэффициентов для отдельных членов суммы. Еще одно свойство метода заключается в том, что для его стабильной работы требуется наличие значительного негравитационного ускорения, хотя это характерно для всех методов коррекции подобного рода.

На практике решаются совершенно аналогичные задачи определения ориентации по измерениям блока определения координат звезд (БОКЗ) на орбите. В противоположность полученному методу, демонстрирующему быстродействие, методы представленные в [4] используют глубокий анализ данных звездных датчиков, гироскопических измерителей угловых скоростей и показаний АСН. При этом для оценки кватерниона ориентации в различных постановках используются МНК, метод множителей Лагранжа, дискретное преобразование Фурье. Зачастую минимизация полученного функционала требует использования итерационного метода Гаусса–Ньютона. Основным недостатком подробного анализа большого объема измерений являются значительные вычислительные затраты, хотя результатом является высокая точность ориентации. Еще одна практическая задача об ускоренном построении орбитальной ориентации по показаниям измерителей угловой скорости и априорном знании об ориентации ракеты-носителя при разделении представлена в [12]. При этом на борту совершается интегрирование назад кинематических и динамических уравнений углового движения, и общее затрачиваемое время на построение ориентации составляет примерно 600 с. Дополнительно результат работы алгоритма можно сравнивать с показаниями прибора ориентации по Земле. Использовать подобные операции во время быстропротекающего процесса спуска длительностью порядка 700 с не представляется возможным, а использование датчиков ориентации по Земле и других оптико-электронных приборов невозможно в силу помех, создаваемых воздушным потоком.

5. Заключение. В работе представлено решение проблемы определения ориентации КА на спуске с помощью МНК. Задача характеризуется высокими требованиями к быстродействию. Хотя для навигационного обеспечения обычно используются традиционные алгоритмы динамической фильтрации, заявленный метод позволил прийти к математически законченной формуле для кватерниона ориентации. Несмотря на то, что исходная СДУ БИНС имеет 10 уравнений и является нелинейной, использование МНК сводит задачу определения ориентации к СЛАУ порядка 3. Это удовлетворительный результат, поскольку решение физической задачи с помощью хорошо изученных математических алгоритмов является важным для систем реального времени. Привести задачу к данной форме позволили исключительные алгебраические свойства кватернионов.