«Применение фильтрации Калмана в задачах определения вращательного движения спутников ...»

Федеральное государственное бюджетное образовательное учреждение

высшего профессионального образования

"Московский государственный технический университет

имени Н. Э. Баумана"

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

Панкратов Владимир Александрович

Применение фильтрации Калмана в задачах

определения вращательного движения спутников

01.02.01 – Теоретическая механика

ДИССЕРТАЦИЯ

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

Научный руководитель д. ф.-м. н., проф., чл.-корр. РАН Крищенко Александр Петрович

Научный консультант д. ф.-м. н., проф.

Сазонов Виктор Васильевич Москва – Содержание

Введение Глава 1. Определение вращательного движения спутника по дан­ ным бортовых измерений вектора напряженности магнитного

поля Земли 1.1. Введение................................ 1.2. Математическая модель вращательного движения спутника, ис­ пользуемая при обработке магнитных измерений......... 1.3. Реконструкция неуправляемого движения методом наименьших квадратов................................ 1.4. Примеры реконструкции неуправляемого движения........ 1.5. Фильтр Калмана............................ 1.6. Примеры реконструкции движения с помощью фильтра Калмана 1.7. Сглаживающий фильтр Калмана.................. 1.8. Примеры сглаживания........................ Глава 2. Определение вращательного движения спутника по дан­........ ным измерений МПЗ и вектора угловой скорости 2.1. Введение................................ 2.2. Измерения угловой скорости на спутнике....... Фотон М- 2.3. Кинематическая модель движения спутника............ 2.4. Методика реконструкции вращательного движения спутника, по данным измерений его угловой скорости и вектора напряженно­ сти МПЗ................................ 2.5. Примеры реконструкции неуправляемого движения........ 2.6. Фильтр Калмана............................ 2.7. Примеры реконструкции неуправляемого движения с помощью фильтра Калмана........................... Глава 3. Проверка согласованности данных измерений магнито­................ метров, установленных на борту ИСЗ 3.1. Введение................................ 3.2. Методика проверки согласованности данных измерений борто­ вых магнитометров.......................... 3.3. Примеры проверки согласованности показаний магнитометров.................................... Заключение Литература Диссертация посвящена задачам реконструкции фактического вращатель­ ного движения искусственных спутников Земли (ИСЗ) научного назначения по данным измерений бортовых датчиков. Основное внимание уделяется ре­ конструкции, необходимой для анализа остаточных микроускорений, которые имели место во время проведения космических экспериментов. Ряд экспери­ ментов по материаловедению, физике жидкости, биологии и медицине весьма чувствительны к остаточным микроускорениям на борту ИСЗ. По этой при­ чине информация о микроускорениях важна для интерпретации получаемых результатов. Анализ многих экспериментов такого рода требует знания толь­ ко квазистатической составляющей микроускорения с частотами менее 0.01 Гц.

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

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

Поскольку многие космические эксперименты с гравитационно-чувстви­ тельными процессами выполняются в течение продолжительного времени, в диссертации предложены методы, позволяющие реконструировать вращатель­ ное движение спутника на интервалах времени до 1 сут. Реконструкция стро­ ится в виде решений динамических и кинематических уравнений движения твердого тела с помощью различных статистических методик. Основное вни­ мание уделено методикам, основанным на фильтрации Калмана. Рассматрива­ ются также интегральные статистические методики, непосредственно исполь­ зующие метод наименьших квадратов. Они используются для проверок, кроме того, некоторые их составные части являются общими с калмановскими мето­ диками. Обычно фильтр Калмана используется для определения движения кос­ мических аппаратов и других механических систем в реальном времени [1, 2].

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

Данные измерений бортовых датчиков, рассматриваемые в диссертации, это — измерения векторов напряженности магнитного поля Земли (МПЗ) и уг­ ловой скорости спутника. Они выполняются с различными шагами по времени, значения которых, как правило, лежат в пределах 5 15 c. Это — косвенные измерения, обработка которых в ряде случаев (например, когда имеются из­ мерения только одного вида) требует применения довольно сложных матема­ тических моделей, основанных на полных (динамических и кинематических) уравнениях движения. По этой причине в диссертации рассмотрены задачи ве­ рификации математических моделей, используемых при обработке косвенных измерений, а также задача проверки показаний бортовых магнитометров.

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

Сначала по измерениям бортовых датчиков, главным образом магнитометров, полученным на некотором отрезке времени, строилась реконструкция враща­ тельного движения спутника на этом отрезке. Затем вдоль найденного дви­ жения микроускорение в заданной точке борта рассчитывалось по известной формуле в функции времени. Реконструкция движения строилась методом наи­ меньших квадратов с использованием решений полных уравнений вращатель­ ного движения спутника. Отрезки времени, на которых эти уравнения позво­ ляли построить адекватную реконструкцию, имели длину от одного до пяти орбитальных витков. Эта длина возрастала вместе с модулем угловой скоро­ сти спутника. Чтобы получить представление о микроускорениях и движении спутника в течение всего полета, движение реконструировалось на нескольких десятках таких отрезков. В данной работе предлагаются несколько вариантов методики реконструкции движения, пригодной для отрезка произвольной дли­ ны. Методика построена на основе фильтра Калмана. Предварительно описан новый вариант методики реконструкции неуправляемого вращательного дви­ жения спутника по магнитным измерениям методом наименьших квадратов, который существенно использован при построении фильтра Калмана. Приво­ дятся результаты сравнения обеих методик на данных, полученных в полете Фотона М-3.

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

В данной главе описаны два подхода к реконструкции вращательного дви­ жения спутника: интегральная статистическая методика, использующая метод наименьших квадратов, и методика, использующая фильтрацию Калмана. В рамках первой методики данные измерений обоих типов, собранные на некото­ ром отрезке времени, обрабатываются совместно. Измерения угловой скорости даны с шагом 12 с на отрезках времени длиной 83 мин. Эти данные сглажива­ ются тригонометрическими полиномами, которые подставляются в кинематиче­ ские уравнения для компонент кватерниона, задающего ориентацию связанной со спутником системы координат относительно гринвичской системы. Получен­ ные таким образом уравнения представляют собой кинематическую модель вра­ щательного движения спутника. Измерения МПЗ выбираются внутри отрезка времени, на котором определены эти уравнения. Их решение, доставляющие ми­ нимум функционалу метода наименьших квадратов считается реконструкцией фактического движения.

Методика, использующая фильтрацию Калмана, разработана в диссерта­ цию на перспективу. В дальнейшем, полет спутников научного назначения, создаваемых ФГУП "ЦСКБ-Прогресс", будет ориентированным. Ориентация спутников солнечными батареями на Солнце будет поддерживаться двигателя­ ми маховиками или гиродинами. В такой ситуации методика мониторинга, ос­ нованная на кинематических уравнениях, станет основной. Для нее возникнет ситуация, упомянутая при описании задачи в первой главе.

Фильтрация организована следующим образом. Заданы две временные сет­ ки. На одной из них заданы измерения угловой скорости. Эта сетка равномерная с шагом 12 с. На второй неравномерной сетке заданы измерения МПЗ. Они зада­ ны с шагом 515 c. Кинематические уравнения интегрируются методом второго порядка точности, специально предназначенным для интегрирования кинемати­ ческих уравнений в кватернионной форме. Интерполяция решения этих уравне­ ний между узлами сетки выполняется с точностью первого порядка и также в кватернионной форме. Интерполяция используется для вычисления расчетных аналогов данных измерений. Обработка магнитных измерений выполняется по одной из схем первой главы.

В рассматривается задача проверки согласованности по­ казаний различных магнитометров, измеряющих МПЗ внутри спутника. Маг­ нитные измерения на начали проводиться не для изучения МПЗ, а с утилитарной целью контроля среды внутри капсулы с научным оборудовани­ ем. Для выполнения измерений использовались несколько трехкомпонентных — 4. Оказалось, что магнитное поле внутри капсулы довольно мало отлича­ ется от МПЗ, и магнитные измерения можно использовать для реконструк­ ции вращательного движения спутника и расчета остаточных микроускорений.

Первые расчеты такого рода были выполнены по измерениям, полученным на лись уже в основном для мониторинга микрогравитационной обстановки. Пере­ численные спутники имели на борту обширные токовые системы, которые вно­ сили заметные возмущения в измерения МПЗ. На каждом из этих трех спут­ ников показания одного или двух магнитометров существенно отличались от расчетных значений МПЗ, показания остальных магнитометров имели значи­ тельные постоянные смещения. Чтобы установить, показания каких магнито­ метров можно использовать для реконструкции вращательного движения спут­ ника, проводились специальные проверки полученных данных. Методика одной из таких проверок описана в диссертации.

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

Приведены примеры применения этой методики при обработке данных, полу­ ченных на Фотоне М-3.

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

Определение вращательного движения спутника по данным бортовых измерений вектора напряженности магнитного поля Земли 1.1. Введение Расчет квазистатических микроускорений на последних четырех спутни­ ках серии (1997 — 2007 гг.) проводился по единой схеме, которая совер­ шенствовалась от экспедиции к экспедиции и описана в [3–14]. Существенную часть этой схемы составляет методика реконструкции фактического вращатель­ ного движения спутника. В методике удачным образом сочетаются два обстоя­ тельства. Во-первых, она основана на полных уравнениях движения ИСЗ. Во­ вторых, квазистатические микроускорения на борту неуправляемого низкоорби­ тального ИСЗ описываются простой формулой, для расчетов по которой надо знать только движение спутника. В каждый момент времени надо знать его радиус-вектор, ориентацию, скорость, угловую скорость и угловое ускорение.

Полные уравнения движения позволяют находить все перечисленные величины.

Кроме того, они позволяют реконструировать фактическое движение спутника по косвенным измерениям. Вращательное движение указанных выше Фотонов было реконструировано по измерениям бортовых магнитометров [4, 7, 11], дат­ чиков угловой скорости [3, 8, 12] и акселерометров [3, 6, 10, 12, 14].

При расчете микроускорений сначала по измерениям бортовых датчиков, полученным на представляющем интерес отрезке времени, строилась рекон­ струкция реального вращательного движения ИСЗ на этом отрезке. Затем вдоль найденного движения вычислялись 13 скалярных функций, позволявших найти квазистатическое микроускорение в любой точке борта в функции време­ ни. Реконструкция движения строилась методом наименьших квадратов. Отрез­ ки времени, на которых используемые уравнения движения позволяли адекват­ ную реконструкцию, имели длину от одного до пяти орбитальных витков. Эта длина возрастала вместе с модулем угловой скорости ИСЗ. Чтобы определить микроускорения в течение всего полета, движение было реконструировано на нескольких десятках таких отрезков, причем соседние отрезки имели перекры­ тие 10 мин. Для спутника полученные результаты описаны в [13].

Реконструкция движения этого спутника в течение 11 суток неуправляемого полета была выполнена по измерениям магнитного поля Земли (МПЗ) и пред­ ставлена 57-ю отрезками. На пересечениях отрезков результаты реконструкции достаточно точно совпали, и с помощью простых вычислительных приемов бы­ ло обеспечено их гладкое сопряжение.

Несколько менее точную методику реконструкции движения, но пригод­ ную для отрезка произвольной длины можно построить на основе фильтра Калмана [15, 16]. Применение этого фильтра для реконструкции вращательно­ го движения спутника по измерениям МПЗ описано в большом числе работ, в основном, зарубежных авторов. В [17] такая реконструкция была использована для расчета квазистатических микроускорений в управляемом движении спут­ ника. Ниже рассматривается применение фильтра Калмана для реконструкции неуправляемого вращательного движения спутников на продолжитель­ ных отрезках времени. Предварительно описан новый вариант методики рекон­ струкции неуправляемого вращательного движения спутника по магнитным из­ мерениям методом наименьших квадратов. Приводятся результаты сравнения обеих методик на данных, полученных в полете Фотона М-3.

Результаты первой главы опубликованы в работе [18].

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

Система 1 2 3 образована главными центральными осями инерции спут­ ника. Точка — его центр масс. При отсутствии специальных указаний компо­ ненты векторов и координаты точек указываются в этой системе.

1 2 3 — гринвичская система координат [19]. Точка — центр Земли, плоскость 1 2 совпадает с плоскостью экватора, положительная полуось пересекает гринвичский меридиан, ось 3 направлена к Северному полюсу.

1 2 3 — квазиинерциальная система координат. Ось 2 направлена вдоль вектора кинетического момента орбитального движения спутника, ось 3 лежит в плоскости 1 2 и направлена в восходящий узел орбиты. Плос­ кость 1 3 совпадает с плоскостью оскулирующей орбиты спутника. Абсо­ лютная величина угловой скорости системы 1 2 3 не превышает нескольких градусов в сутки.

Матрицы перехода от системы 1 2 3 к системам 1 2 3 и 1 обозначим соответственно,=1 и,=1. Здесь и — косинусы углов, которые ось образует с осями и. Примем следующие способы параметризации этих матриц.

Положение системы 1 2 3 относительно гринвичской системы коор­ динат будем задавать нормированным кватернионом Q = (0, 1, 2, 3 ), 2 + 2 + 2 + 2 = 1. Элементы матрицы выражаются через ком­ поненты Q с помощью формул [20, 21] Кватернионный вид формул перехода Положение системы 1 2 3 относительно системы 1 2 3 будем зада­ вать углами, и, которые введем следующим образом. Если совместить точки и, то система 1 2 3 может быть переведена в систему 1 тремя последовательными поворотами: 1) на угол + /2 вокруг оси 2, 2) на угол вокруг новой оси 3, 3) на угол вокруг новой оси 1, совпадающей с осью 1. Элементы матрицы выражаются через эти углы с помощью формул Углы, и используются для графического представления вращатель­ ного движения спутника — это движение удобно иллюстрировать графиками зависимости указанных углов от времени. Кватернион Q входит в фазовый вектор вращательного движения спутника.

Уравнения движения спутника состоят из двух подсистем. Одна подси­ стема описывает движение центра масс спутника, другая — его вращательное движение. Подсистема уравнений движения центра масс записывается в грин­ вичской системе координат [6, 7]. В ней учитываются нецентральность гравита­ ционного поля Земли и сопротивление атмосферы. Нецентральность поля учи­ тывается с точностью до членов порядка (16,16) включительно в разложении гравитационного потенциала Земли в ряд по шаровым функциям. Атмосфера считается вращающейся вместе с Землей, ее плотность рассчитывается согласно модели ГОСТ Р 25645.166-2004. Параметры атмосферы и баллистический коэф­ фициент спутника считаются неизменными на всем интервале интегрирования уравнений движения.

Подсистема уравнений вращательного движения образована динамически­ ми уравнениями Эйлера для компонент угловой скорости спутника и кинема­ тическими уравнениями для компонент кватерниона Q [20, 21]. В уравнениях Эйлера учитываются гравитационный и восстанавливающий аэродинамический моменты, а также учитывается гиростатический момент внутренних устройств спутника (вентиляторов, роторов и т. п.). Подсистема уравнений вращательного движения имеет вид Здесь, и — компоненты вектора абсолютной угловой скорости спутника, геоцентрического радиуса-вектора точки и скорости этой точки относительно гринвичской системы координат, — моменты инерции спутника относитель­ но осей, — параметры аэродинамического момента, — отнесенные к 1 компоненты гиростатического момента внутренних устройств спутника, — угловая скорость вращения Земли, — масштабирующий множитель. При численном интегрировании уравнений (1.1) единицами измерения времени и длины служат 103 с и 103 км, единицы измерения других величин: [ ] = км/с, [ ] = [ ] = 103 1, [ ] = см/кг, [ ] =кг/м3, = 1010.

В дальнейшем используется более компактная запись кинематических уравнений [22] где — символ Леви-Чивиты (равен 1, если,, — четная перестановка чисел 1, 2 и 3, равен 1 для нечетной перестановки и равен 0 в остальных случаях). Переменные зависимы — связаны условием нормировки кватерни­ она Q. Если это условие выполнено в начальный момент времени, то в силу свойств уравнений (1.2) оно будет выполняться тождественно. Следовательно, достаточно обеспечить условие нормировки только в начальный момент.

Параметры и в уравнениях (1.1) считаются известными. Для Фотона М- на каждом интервале обработки данных измерений (см. ниже), но их значения определяются в результате этой обработки наряду с неизвестными начальными условиями движения спутника, т. е. и служат параметрами согласования.

1.3. Реконструкция неуправляемого движения методом наименьших квадратов Методику реконструкции неуправляемого вращательного движения низко­ орбитального КА по магнитным измерениям опишем на примере Фотона M-3.

На борту этого спутника находились четыре трехкомпонентных магнитометра, входящих в состав аппаратуры DIMAC [11]. Аппаратура предназначалась для измерения микроускорений на борту спутника. Основными ее датчиками бы­ ли акселерометры. Магнитные измерения проводились для реконструкции вра­ щательного движения спутника с целью проверки показаний низкочастотного акселерометров расчетным путем [9, 10].

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

Измерения разных магнитометров оцифровывались на единые моменты време­ ни, промежутки между которыми варьируются в пределах от 1 до 12 с, а в среднем составляют около 5 с. Для реконструкции движения брались сплош­ ные ряды этих измерений, охватывающие интервалы времени длиной от 2 до 8 ч. Выбранные данные представляли собой совокупность чисел где ( = 1, 2, 3) измеренные значения компонент векторана пряженности магнитного поля в момент, < < · · · <. Полагаем, что эти компонен­ ты с точностью до постоянных смещений, а также малых ошибок измерений и координатных преобразований совпадают с компонентами напряженности МПЗ в системе координат 1 2 3.

Следуя методу наименьших квадратов, аппроксимацией фактического дви­ жения спутника на отрезке будем считать решение системы (1.1), доставляющее минимум функционалу Здесь — оценки постоянных смещений в измерениях, () — расчетные значения компонент напряженности МПЗ в гринвичской системе координат в момент времени. Функции () строятся вдоль известной орбиты спутника с использованием аналитической модели МПЗ IGRF2005.

Функционал (1.4) получен в результате преобразования стандартного функционала метода наименьших квадратов, возникающего при уравнивании нимизация проводится по начальным условиям решения ( ), ( ) и параметрам математической модели, ( = 1, 2, 3; = 0, 1, 2, 3). При этом учитывается условие нормировки Для простоты письма все уточняемые величины объединим в один вектор R13. В принятых обозначениях = (), * = argmin () — искомая оценка вектора. Минимизация () выполнялась в несколько этапов разными методами. Ее описание начнем с заключительного этапа, на котором применял­ ся метод Гаусса — Ньютона [23].

На каждой итерации этого метода поправки ( ) к имеющимся значе­ ния ( ) ищутся в виде (ср. уравнения (1.2)) Параметры суть компоненты вектора бесконечно малого поворота [24], зада­ ющего изменение ориентации спутника в окрестности положения Q( ). Эти параметры и поправки ( ),, находятся из системы нормальных уравнений с матрицей 12 и правой частью 12 :

Здесь 1, 2,..., 12 — обозначения величин 1, 2, 3, 1 ( ), 2 ( ), 3 ( ), 1, 2, 3, 1, 2, 3 в указанном порядке, ()/ — псевдопроизводные, служащие для представления истинных производных Псевдопроизводная — это не частная производная некоторой функции по ка­ кому-то параметру. Запись ее в виде частной производной используется лишь для удобства. Такую запись следует воспринимать как единый символ с двумя индексами. В кинематике твердого тела угловая скорость служит для расче­ та производных по времени, а псевдопроизводная — для расчета производных по параметру (ср. выписанные выражения для / c уравнениями (1.2) и формулами (1.6)). В обозначении / индекс указывает векторную ком­ поненту, индекс — номер параметра, по которому выполняется дифференциро­ вание. Дифференцируя уравнения (1.7) и уравнения вращательного движения твердого тела, записанные в кватернионной форме, получим С учетом уравнений вращательного движения твердого тела, записанных в ква­ тернионной форме, последнее равенство примет вид Это равенство обеспечивается при равенстве нулю одного из сомножителей. Из равенства Q 1, следует запись системы дифференциальных уравнений для определения значений псевдопроизводных в координатном виде Эти уравнения интегрировались совместно с уравнениями (1.1) и уравнениями в вариациях для /. Последние получаются дифференцированием по первых трех уравнений (1.1), причем производные / выражаются через / с помощью приведенных выше формул. Ненулевые начальные условия для / и / имеют вид Прибавление найденных поправок ( ) к имеющимся значениям ( ) на­ рушает условие (1.5), поэтому новый кватернион ориентации нормируется. Вне­ сенные нормировкой изменения уточненных компонент кватерниона являются величинами второго порядка относительно ( ).

Интегрирование уравнений (1.1) и указанных выше уравнений в вариаци­ ях выполняется одним из методов Дормана–Принса 8-го порядка. Это — метод типа Рунге–Кутты, реализованный в стандартной процедуре DOP853 [25]. Ме­ тод и программа позволяют строить полином, интерполирующий вычисляемое решение внутри шага интегрирования. Этот полином используется для вычисле­ ния функционала (1.4), матрицы и правой части системы нормаль­ ных уравнений. При этом интегрирование уравнений движения и уравнений в вариациях выполняется с оптимальным достаточно большим шагом, величина которого выбирается по критерию локальной точности интегрирования. Нали­ чие интерполяционного полинома и возросшее быстродействие персональных компьютеров позволили осуществлять совместную обработку всех собранных данных измерений, не проводя их предварительную обработку. Это — одно из отличий описываемой методики от методик, использованных в [4, 7, 11]. Дру­ гое отличие заключается в использовании компонент кватерниона Q в качестве кинематических переменных уравнений движения спутника и в способе уточне­ ния начальных условий этих переменных при реализации метода Гаусса–Нью­ тона. В [4, 7, 11] кинематическими переменными служили величины 1, ( = 1, 2, 3), начальные условия которых параметризовались тремя углами точ­ но так же, как величины 1, 2 параметризуются углами, и.

Вернемся к описанию методики. Точность аппроксимации измерений и оценки * будем характеризовать, следуя методу наименьших квадратов, со­ ответствующими стандартными отклонениями [26]. Последние вычисляются в предположении, что ошибки в измерениях некоррелированы и имеют одинако­ вые дисперсии, средние значения ошибок в измерениях одной и той же вектор­ ной компоненты напряженности МПЗ равны [27]. Такой подход выбран из сооб­ ражений удобства и вида функционала (1.4). При сделанных допущениях * — случайный вектор, который имеет приблизительно нормальное распределение со средним значением, равным истинному значению. Вследствие условия (1.5) это распределение — несобственное, т. е. имеет вырожденную ковариационную матрицу. Чтобы избежать вырождения и сделать характеризацию ошибок более наглядной, ошибки ( ) в задании компонент ( ) вектора * предста­ вим в виде (1.6), где теперь образуют случайный вектор бесконечно малого поворота. Величины имеют нулевые математические ожидания и вместе с ошибками остальных компонент * описываются ковариационной матрицей Здесь 2 оценка дисперсии ошибок в измерениях, — матрица, вычис­ ленная в точке *. Точность аппроксимации измерений будем характеризовать стандартным отклонением, точность оценки * — стандартными отклонения­ ми ( = 1, 2,..., 12). Стандартные отклонения величин, (0 ), и будем обозначать,,,.

Чтобы обеспечить надежную сходимость описанного процесса, надо преду­ смотреть возможность его регуляризации и иметь достаточно точное началь­ ное приближение точки *. Регуляризация сводилась к предварительному ис­ пользованию метода Левенберга–Марквардта [23] перед переходом к методу Гаусса–Ньютона [28]. В большинстве случаев метод Левенберга–Марквардта плавно трансформировался в метод Гаусса–Ньютона, но иногда после оконча­ ния его работы метод Гаусса–Ньютона расходился. В таком случае в качестве * принимался результат, полученный методом Левенберга–Марквардта. Резуль­ тат полученный таким способом обеспечивал требуемую точность реконструк­ ции движения спутника.

Поиск начального приближения точки * выполнялся так. Перед миними­ зацией функционала (1.4) по 12 параметрам, выполнялась его минимизация по 9 параметрам при 1 = 2 = 3 = 0. Она выполнялась сначала методом слу­ чайного поиска, затем методом Левенберга–Марквардта и, наконец, методом Гаусса–Ньютона. В ряде случаев даже применение такой редуцированной моде­ ли движения — системы (1.1) при 1 = 2 = 3 = 0 — обеспечивало достаточно высокую точность реконструкции движения спутника.

Начальная точка для минимизации функционала (1.4) по 9 параметрам находилась посредством минимизации так называемого укороченного функци­ онала, который задавался теми же формулами (1.4), но при уменьшенном в несколько раз значении и специальном выборе начальных условий движе­ ния и параметров уравнений (1.1). Сокращенный отрезок данных имел длину 15 — 30 мин. В уравнениях движения принималось = = 0 ( = 1, 2, 3).

Начальные условия для кватерниона задавались в виде Здесь — параметр. Начальные условия для псевдопроизводных по параметру Выписанные соотношения получены при следующих предположениях. Орт e, имеющий в системах координат 1 2 3 и 1 2 3 компоненты и соответ­ ственно, в этих системах фиксирован. В системе 1 2 3 он направлен по рас­ четному вектору напряженности МПЗ, а в системе 1 2 3 — по измеренному вектору. Поворот спутника вокруг этого орта может быть произвольным и зада­ ется углом. Укороченный функционал рассматривался в функции четырех па­ раметров: и ( ) ( = 1, 2, 3). Для его минимизации последовательно приме­ няются случайный перебор в параллелепипеде {0 2, ( ) } ( и заданы), случайный поиск и метод Левенберга–Марквардта.

1.4. Примеры реконструкции неуправляемого движения тодом наименьших квадратов приведены на на рис. 1.1-1.4 (более подробные результаты, полученные с помощью той же математической модели, но другим программным обеспечением представлены в [11, 13]). Реконструкция выполнена на двух интервалах времени. Интервал на рис. 1.1, 1.2 приходится на середину полета, интервал на рис. 1.3, 1.4 — на его конец. Рис. 1.1 и 1.3 иллюстриру­ ют движение спутника относительно квазиинерциальной системы координат 1 2 3, рис. 1.2 и 1.4 характеризуют точность аппроксимации данных изме­ рений. Графики на рисунках построены на отрезках, в подписях к рисункам указаны начальные точки интервалов, значения и оценки стандартного отклонения ошибок в измерениях.

Каждая страница рисунков естественным образом разбивается на две ча­ сти — левую и правую. На рис. 1.1 и 1.3 в левых частях приведены гра­ фики зависимости от времени углов, и, а также график разности () = ()0 1 ( ), где 0 +1 ( ) — линейная аппроксимация функ­ ции (), построенная методом наименьших квадратов. В правых частях этих рисунков помещены графики компонент угловой скорости () в найденных решениях уравнений (1.1).

В левых частях рис. 1.2 и 1.4 приведены графики функций () (см. (1.4)) и ломаные, проходящие через точки (, ломаная и график аппроксимирующей ее функции () изображены в единой системе координат. Ломаные и аппроксимирующие их графики практически совпадают, поэтому в правых частях рис. 1.2 и 1.4 приведены ломаные, прохо­ Приведем некоторые числовые характеристики найденных реконструкций.

В решении на рис. 1.1 и 1. Стандартные отклонения оценок начальных условий и уточняемых параметров составляют В решении на рис. 1.3 и 1. Стандартные отклонения здесь выражены в радианах, стандартные отклоне­ ния остальных величин — в единицах, в которых интегрируюся уравнения (1.1), в частности, [ ] = 103 1.

Оценки параметров в решении на рис. 1.3 и 1.4 выглядят неправдоподоб­ но большими. Это связано с тем, что рассматриваемое решение построено для продолжительного отрезка времени, на котором начинает сказываться неадек­ ватность модели. Если длину отрезка уменьшить, то оценки параметров уменьшатся.

Рис. 1.1. Движение реконструированное методом наименьших квадратов. Момент = 0 на графиках соответствует 22:46:42.00 ДМВ 17.09.2007.

Рис. 1.2. Аппроксимация магнитных измерений, полученная методом наименьших квадратов, = 2340, = 958. Момент = 0 на графиках соответствует 22:46:42.00 ДМВ 17.09.2007.

Рис. 1.3. Движение Фотона реконструированное методом наименьших квадратов. Момент = 0 на графиках соответствует 09:09: ДМВ 22.09.2007.

Рис. 1.4. Аппроксимация магнитных измерений, полученная методом наименьших квадратов, = 4000, = 831. Момент на графиках соответствует 09:09:05 ДМВ 22.09.2007.

Чтобы получить более содержательную оценку точности определения дви­ жения спутника с помощью описанной методики, была проведена обработка измерений на 11 пересекающихся интервалах времени, охватывающих 16.7 ч.

Каждый интервал имел длину примерно 6200 с. Смежные интервалы имели пересечение около 800 с. Начальные точки интервалов приведены в табл. 1.1.

Реконструкция движения спутника на всех этих интервалах (11 решений урав­ нений (1.1)) представлена на рис. 1.5. В левой части рисунка приведены гра­ фики зависимости от времени углов, и, а также график разности () = () 0 1 ( ), где 0 + 1 ( ) — линейная аппроксима­ ция функции (), построенная методом наименьших квадратов. В средней и правой частях рисунка приведены графики зависимости от времени компонент угловой скорости и кватерниона ориентации. Как видно из рисунка, на стыках интервалов разные решения несколько отличаются, но в целом реконструкция выглядит плавной. Оценки параметров модели на этих интервалах приведены в табл. 1.1.

В [13] приведена аналогичная реконструкция всего неуправляемого поле­ та Там использовались более сложные уравнения движения — дополнительно учитывался действующий на спутник магнитный момент, пара­ метры которого уточнялись, и стыковка смежных решений получилась более точной. Следует отметить, что и интервалы времени, на которых выполнялась реконструкция в [13] были длиннее — в основном 2 3 орбитальных витка.

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

Рис. 1.5. Результаты реконструкции движения на 11 смежных промежутках времени.

Таблица 1.1. Параметры модели, основанной на динамических уравнениях вращательного движения 10 0.7334 0.0067 0.2697 0.0030 0.1118 0. 11 0.8436 0.0070 0.2702 0.0041 0.0910 0. 1.5. Фильтр Калмана Дискретный фильтр Калмана в достаточно общем виде можно описать следующим образом [15, 16, 25]. Рассмотрим механическую систему, некоторые параметры движения которой измеряются в дискретные моменты времени. Ма­ тематическую модель движения системы и процесса измерений примем линей­ ной и представим в виде Здесь — фазовый вектор системы, — вектор измерений, — вектор случайных ошибок измерений, — вектор случайных возмущений системы.

Размерность векторов, и порядок квадратных матриц считаем не за­ висящими от, размерности и размеры остальных векторов матриц полагаем согласованными должным образом. Полагаем также, что все матрицы имеют полный ранг. В частности, все матрицы — невырожденные. Полагаем далее, что векторы и некоррелированы при любых значениях и, векторы и, а также и, некоррелированы при любых =. Математические ожидания и ковариационные матрицы этих векторов определяются соотноше­ ниями где и — симметричные, положительно определенные матрицы. Задача состоит в том, чтобы по данным измерений 1, 2,..., найти оценку вектора и ковариационную матрицу ошибки.

Предположим сначала, что введенные выше случайные величины и имеют нормальное распределение. Тогда фазовый вектор и его оценка также будут нормальными случайными величинами и решение поставленной задачи об оценивании движения системы можно свести к задаче нахождения условной плотности распределения ( | ). Здесь = {1, 2,..., }. В этом Вывод рекуррентных соотношений для расчета величин и можно осуществить посредством формальных преобразований условных плотностей вероятностей [29]. С учетом равенства (, ) (, 1, ) имеем ( |, 1 ) = ( | ), справедливого в силу второго соотношения (1.9).

В силу первого соотношения (1.9) ( |1, 1 ) = ( |1 ), поэтому и, следовательно, Здесь интегрирование выполняется по всему фазовому пространству системы.

В силу нормальности рассматриваемых случайных величин [30] (здесь и ниже несущественные множители опускаем) Посредством несложных, но довольно громоздких преобразований подынте­ гральное выражение в последнем интеграле можно преобразовать к виду (чле­ ны, не содержащие и 1, опущены) Вычислив выписанный интеграл, получим Согласно второму соотношению (1.9) и принятым допущениям Перемножив правые части двух последних формул, получим после преобразо­ ваний Здесь Соотношения (1.10) задают рекуррентные формулы для расчета величин и [31, 32]. Отыскание величин 0 и 0 — отдельная задача. Полагаем, что она решена. В ряде случаев величины 0 и 0 можно выбирать достаточно произвольно; с увеличением рекуррентные формулы (1.10) начинают давать правильные значения и.

Как следует из вывода формул (1.10), вектор = 1 — несмещенная оценка вектора по измерениям 1, — ковариационная матрица этой оценки.

Соотношения (1.10) позволяют находить оценку фазового вектора систе­ мы по мере поступления измерений. Обычно их применяют в случае, когда dim < dim. Если выполнено обратное неравенство и dim зависит от, то вследствие необходимости обращать матрицу порядка dim форму­ лы (1.10) следует преобразовать к более удобному виду. На промежуточном этапе вывода формул (1.10) возникают соотношения [33, 34] Фильтр, построенный на их основе, оперирует с матрицами 1, 1, 1,, 1. Обращаются только матрицы (для вычисления ), и порядка dim. Матрицы, 1 являются исходными вместо,.

В задачах обработки магнитных измерений на спутниках Земли матрица близка к единичной, а элементы существенно меньше 1. С учетом Рассмотрим последнюю формулу в случае, когда =, — поло­ жительный скаляр, — единичная матрица. Такое предположение означает, в частности, что единицы измерения компонент фазового вектора выбраны под­ ходящим образом). Тогда 1 =. Если и — собственные вектор и число матрицы ( = 0), то и 2 — собственные вектор и число матрицы 1. Матрицы 1 и имеют одинаковые собственные векторы, а их собственные числа в случае 1 близки. В такой ситуации можно взять 1 = при подходящем (0, 1), 1.

Выписанные выше соотношения относятся к линейному случаю. Перейдем теперь к нелинейной системе. Математическую модель ее движения и процесса измерений примем в виде В выписанных соотношениях () и () — гладкие функции, матрицы Яко­ би которых имеют полный ранг; векторы,, и имеют прежний смысл, причем ошибки и имеют указанные выше первые и вторые моменты.

Рассмотрим получение оценки вектора по измерениям 1, 2,...,.

Эту оценку и ее ковариационную матрицу по-прежнему будем обозначать и. Рекуррентные соотношения, позволяющие находить искомые величины по­ лучаются эвристической модификацией соотношений (1.11). А именно, оценка в (1.11) минимизирует по выражение Иными словами, эта оценка является оценкой метода наименьших квадратов.

Второе слагаемое в выписанном выражении отвечает за ошибки измерений, пер­ вое слагаемое учитывает априорную информацию о значении. По аналогии и с учетом малости случайных величин в модели (1.12) оценку вектора в рамках этой модели будем искать из условия минимума функции где теперь Описанный переход от линейного фильтра к нелинейному не единствен.

Приведем пример другого перехода. Минимизация по функции 1 ( ) экви­ валентна минимизации по 1 функцию Первое слагаемое в правой части формулы для 1 [ (1 )], отвечающее за априорную информацию с высокой точностью можно представить в виде Здесь матрица имеет тот же вид, что и в формуле для 1 ( ). Положим и примем в качестве оценки величину = [arg min 2 (1 )]. Понятно, что эта оценка отличается от предыдущей оценки = arg min ( ), но при определенных условиях они будут близки.

Минимизацию функций 1 ( ) и 2 (1 ) можно выполнять методом Гаусса–Ньютона, в качестве начальных приближений использовать точки и 1 соответственно. При минимизации функции 1 ( ) нормальные уравнения имеют вид Здесь * — имеющееся приближенное значение величины, которму после ^ решения системы (1.13) присваивается значение * +. Итерационное уточ­ нение * прекращается, когда норма поправки станет приемлемо малой.

Тогда принимается = *.

При минимизации функции 2 (1 ) нормальные уравнения имеют вид — функция, обратная функции (). Итерационное уточнение * выполня­ ется аналогично только что описанному уточнению оценки. По окончании итераций принимается = (* ).

Нелинейные фильтры, основанные на решении систем (1.13), (1.14) бы­ ли использованы для обработки магнитных измерений, выполненных на В программных реализациях этих фильтров под уравнениями дви­ Фотоне М-3.

жения спутника понимается совокупность уравнений (1.1) и уравнений Для ссылок такую систему обозначим (1.1 ). Компонентами ее фазового векто­ ра R16 являются величины,,, ( = 1, 2, 3) и ( = 0, 1, 2, 3).

Значения фазового вектора суть значения переменных уравнений движения спутника в узлах временной сетки { } ( = 0, 1, 2,... ), 0 < 1 < 2 <...

В реализации, использующей минимизацию функции 2 (1 ), эта сетка равно­ мерная: = 0 +, = 200 500с. Выбираемое значение в несколько раз превышает длину оптимального шага интегрирования уравнений (1.1) с помо­ щью процедуры DOP853. В реализации, основанной на минимизации функции 1 ( ), узлы являются концевыми точками этих оптимальных шагов. Здесь +1 = 50 100с.

В обеих реализациях формирование векторов измерений и вычисле­ ние функций (), () выполняется следующим образом. Решение систе­ мы (1.1 ) с начальным условием ( ) = обозначим = (,, ). Тогда () = (, 1, ). Измерения напряженности магнитного поля (ср. (1.3)), (1, ] — объединим в вектор. Соответствующие расчетные аналоги этих измерений + ( ), вычисленные вдоль решения (, 1, ) с использо­ ванием интерполяционного полинома, обозначим (). Имеют место формулы 1 = 1 *, * в системе (1.14) фактически являются эле­ ментами R15 — четыре кватернионные компоненты этих векторов (они малы) заменены тремя компонентами векторов бесконечно малых поворотов. Такая замена была использована в п. 1.3 при формировании системы нормальных уравнений метода наименьших квадратов. По этой причине матрицы, и т.п. имеют порядок 15. Вместо ковариаций компонент кватернионов использу­ ются ковариации компонент векторов бесконечно малых поворотов.

Приведем некоторые детали формирования и решения системы (1.13).

Пусть для точки 1 известны оценка 1 и матрица 1. Сначала посред­ ством интегрирования уравнений (1.1 ) от точки 1 к точке (фактически здесь и ниже интегрируются уравнения (1.1) и описанная в п. 1.3 система урав­ нений в вариациях) находятся прогноз и матрица 1. Расчет этой матри­ цы выполняется по третьей и последующим формулам (1.11), в которых теперь = (^1 )/. Затем реализуется итерационный процесс получения оцен­ ки. Матрица и правая часть системы (1.13) вычисляются посредством инте­ грирования уравнений (1.1 ) от точки к точке 1. При этом часть элементов матрицы и компонент вектора [ ( (1 )] рассчиты­ ваются по формулам, используемым при формировании системы нормальных уравнений метода наименьших квадратов. Матрица 1 на всех итерациях остается неизменной. Ковариационная матрица найденной оценки равна об­ ратной матрице системы (1.13) на последней итерации уточнения * (ср. вторую формулу (1.11), в которой вычисляется при = ). Фактически в фильтре используется только эта матрица системы (1.13), т. е. матрица 1, а матрица служит только для расчета стандартных отклонений оцениваемых величин.

Процесс решения системы (1.14) организован несколько иначе. На каж­ дой итерации уравнения (1.1 ) интегрируются от точки 1 к точке. На пер­ вой итерации вычисляется матрица, которая при последующих итерациях остается неизменной. Часть элементов матрицы и компонент вектора [ (* )] рассчитываются так же, как при формировании системы нормальных уравнений метода наименьших квадратов в п. 1.3. Ковариационная матрица оценки = (* ), * = arg min 2 (1 ) определяется формулой Фактически в фильтре используется матрица Оценка 0 и матрица 01, относящиеся к начальной точке 0, находились методом наименьших квадратов посредством минимизации функционала (1.4) 01 = 2. Однако соотношения, связывающие все рассматриваемые ковари­ ационные матрицы, являются однородными первой степени. Кроме того, рекон­ струируемое движение спутника на довольно продолжительных отрезках сетки { } можно считать установившимся. В такой ситуации удобно в качестве принять единичную матрицу, положить 01 =, взять диагональной мат­ рицей, не зависящей от, а ковариационную матрицу оценки вычислять как 0, где 0 — характерное значение стандартного отклонения ошибок в измерениях отдельных компонент МПЗ в строительной системе координат. В расчетах следующего пункта диагональные компоненты, отвечающие фазо­ вым переменным, составляют 2 · 104, отвечающие параметрам — 2 · 102 (в едини­ цах измерения переменных и параметров системы (1.1)), в качестве 0 берется среднеквадратичная ошибка аппроксимации данных измерений. Эта ошибка и матрицы 0 рассчитываются после реконструкции движения на всем выбран­ ном отрезке времени. Матрицы, оценки и еще ряд величин, относящихся к интервалам (1, ] хранятся в памяти компьютера.

Точность аппроксимации измерений вектора МПЗ на каждом полуинтер­ вале (1, ] характеризуется величинами где — количество элементов множества ().

1.6. Примеры реконструкции движения с помощью фильтра Калмана времени, пересекающемся с объединением 11 интервалов из п. 1.4 (см. табл. 1.1).

Реконструкции, полученные посредством решения системы (1.14) приведены на рис. 1.6–1.7. Рис. 1.6–1.9 получены при +1 = 200с. Левая часть рис. 1. содержит графики зависимости от времени углов (ср. рис. 1.1, 1.3), и, а также график разности () = ()0 1 (0 ), где 0 +1 (0 ) — линейная аппроксимация функции (), построенная методом наименьших квадратов. В средней части рисунка помещены графики компонент угловой скорости (), в левых частях — ломаные, характеризующие ошибки аппроксимации данных измерений (ср. правые графики на рис. 1.2, 1.4). На рис. 1.7 изображены гра­ фики величин,, и. Эти величины — кусочно постоянные функции, на интервалах (1, ] они сохраняют свои значения. На рис. 1.8 и 1.9 приве­ дены графики стандартных отклонений оценок. Каждая компонента оцен­ ки представлена соответствующим стандартным отклонением, за исключением компонент кватерниона. Последние представлены стандартными отклонениями компонент вектора бесконечно малого поворота ( = 1, 2, 3). Это также кусочно постоянные функции с промежутками постоянства (1, ]. В данном примере = 0 = 482.

Реконструкция движения на том же отрезке времени посредством решения системы (1.13) приведена на рис. 1.10–1.13. Эта серия рисунков организована аналогично рис. 1.6–1.9. В данном случае 0 = 1009. Точность аппроксимации измерений получилась несколько хуже, чем при использовании системы (1.14).

Это видно и из сравнения сравнения рис. 1.12, 1.13 с рис. 1.8, 1.9. Указанный эффект обусловлен уменьшением числа измерений на полуинтервале (1, ] при переходе от варианта с системой (1.14) к варианту с системой (1.13) — длина этого интервала стала в 2 3 раза короче. Вариант с системой (1.13) оказался чувствителен к этому числу. Даже кратковременный сбой в поступле­ нии измерений может привести к некоторому нарушению работы фильтра. На рис. 1.8 заметен выброс при 830 мин. Он обусловлен именно малым числом измерений на соответствующем полуинтервале (1, ].

Варьируя элементы матрицы, можно заметно повысить точность ап­ проксимации измерений даже при использовании варианта с системой (1.13).

В частности, в данном примере можно достичь значения 0 < 300. Но такое уменьшение сопровождается увеличением диапазона изменения параметров, в (1.1). Аналогичный результат получается, если вместо использования принять 1 =, 1, < 1 (см. выше). Здесь получается высокая точ­ ность аппроксимации и неоправданно большие (по сравнению с оценками мето­ да наименьших квадратов) значения оцениваемых параметров уравнений (1.1).

Рис. 1.6. Реконструкция движения с помощью фильтра Калмана. Шаг по времени 200с, 0 = 481. Момент = 0 на графиках соответ­ ствует 00:57:00 ДМВ 18.09.2007.

Рис. 1.7. Реконструкция движения с помощью фильтра Калмана. Шаг по времени 200с. Момент = 0 на графиках соответствует 00:57:00 ДМВ 18.09.2007.

Рис. 1.8. Реконструкция движения с помощью фильтра Калмана. Шаг по времени 200с. Момент = 0 на графиках соответствует 00:57:00 ДМВ 18.09.2007.

Рис. 1.9. Реконструкция движения с помощью фильтра Калмана. Шаг по времени 200с. Момент = 0 на графиках соответствует 00:57:00 ДМВ 18.09.2007.

Рис. 1.10. Реконструкция движения с помощью фильтра Калмана. Автоматический выбор шага по времени. Момент = 0 на графиках соответствует 00:57:00 ДМВ 18.09.2007.

Рис. 1.11. Реконструкция движения с помощью фильтра Калмана. Автоматический выбор шага по времени. Момент = 0 на графиках соответствует 00:57:00 ДМВ 18.09.2007.

Рис. 1.12. Реконструкция движения с помощью фильтра Калмана. Автоматический выбор шага по времени. Момент = 0 на графиках соответствует 00:57:00 ДМВ 18.09.2007.

Рис. 1.13. Реконструкция движения с помощью фильтра Калмана. Автоматический выбор шага по времени. Момент = 0 на графиках соответствует 00:57:00 ДМВ 18.09.2007.

1.7. Сглаживающий фильтр Калмана Описанные выше процессы фильтрации позволяют реконструировать вра­ щательное движение спутника по мере поступления данных измерений. Их мож­ но использовать для реконструкции в реальном времени. Для апостериорной реконструкции движения по всей полученной совокупности данных обычно ис­ пользуют другие алгоритмы. Задача апостериорной реконструкции формулиру­ ется следующим образом. Пусть на некотором отрезке времени были получены измерения при = 1, 2,...,, и по мере поступления этих измерений по формулам (1.10) рассчитаны величины,, и. После того как процесс получения измерений был завершен возникает задача уточнения оценок и с использованием всей имеющейся измерительной информации. Уточненную оценку вектора обозначим, ковариационную матрицу новой оценки обо­ значим. Ниже эта задача решается с помощью алгоритма, основанного на результатах [35]. Приведем, следуя [29], вывод основных соотношений этого ал­ горитма в терминах модели (1.9) и предположений п. 1.5. Будем основываться на формальных преобразованиях условных плотностей вероятностей.

Найдем ( | ). По определению плотности условной вероятности имеем Отсюда Далее Здесь в последнем выражении выписаны только члены второй степени по ком­ понентам векторов, +1 и члены первой степени по компонентам вектора. Их знание достаточно, чтобы вычислить интеграл в формуле где интегрирование выполняется по всему фазовому пространству системы.

Рассмотрим вспомогательный интеграл где и постоянные матрицы и векторы. Квадратичную форму (1, 2 ) представим в виде С использованием последней формулы получим Приняв 1 = +1, 2 = и сравнив ln (+1, | ) с (1, 2 ), найдем Согласно установленному соответствию где = 11 12 22 12. Прямые вычисления дают Отсюда При использовании формул (1.17) принимается =, =.

Перейдем к нелинейной модели (1.12). Рекуррентные соотношения, позво­ ляющие осуществлять сглаживание, получаются эвристической модификацией соотношений (1.17). Оценка, выражаемая первой формулой (1.17), миними­ зирует по выражение т.е. является оценкой метода наименьших квадратов. По аналогии и с учетом малости случайных величин в модели (1.9) оценку вектора в рамках этой модели будем искать из условия минимума функции Минимизация функции ( ) выполняется методом Гаусса–Ньютона, началь­ ным приближением служат точка. Нормальные уравнения имеют вид Здесь * — имеющееся приближенное значение величины, которому после решения системы (1.18) присваивается значение * +. Итерационное уточ­ нение * прекращается, когда норма поправки станет приемлемо малой.

Тогда принимается = *.

Процесс вычисления оценок, идет в направлении от точки к точке 1. При переходе от точки +1 к точке матрица и правая часть системы (1.18) находятся интегрированием уравнений (1.1 ) в обратном направлении — от точ­ ки к точке +1. Расчет матрицы 1 выполняется по второй формуле (1.17), где +1 имеет вид (1.18). В начале процесса =, =.

Чтобы реализовать этот процесс, приходится хранить в памяти компью­ тера все вычисленные в процессе фильтрации "вперед"величины, и.

В [35] приведены реккуррентные формулы, которые позволяют вычислить мат­ рицы и в процессе фильтрации "назад". Однако применение этих фор­ мул часто дает отрицательно определенные матрицы (возможно сказывается нелинейность задачи), поэтому их использование непрактично.

1.8. Примеры сглаживания Результаты сглаживания приведены на рис. 1.14–1.21. Серия рис. 1.14–1. получена для варианта, в котором фильтрация "вперед"выполнена с использо­ ванием системы (1.14), серия рис. 1.18–1.21 относится к варианту, в котором фильтрация "вперед" выполнена с использованием системы (1.13). Как пока­ зывает сравнение этих серий рисунков, один и тот же метод сглаживания при использовании разной (хотя очень похожей) исходной информации дает заметно различающиеся результаты. Однако микроускорения для обеих реконструкций оказываются практически одинаковыми.

Рис. 1.14. Сглаживание с помощью фильтра Калмана. Шаг по времени 200с. Момент = 0 на графиках соответствует 00:57:00 ДМВ 18.09.2007.

Рис. 1.15. Сглаживание с помощью фильтра Калмана. Шаг по времени 200с. Момент = 0 на графиках соответствует 00:57:00 ДМВ 18.09.2007.

Рис. 1.16. Сглаживание с помощью фильтра Калмана. Шаг по времени 200с. Момент = 0 на графиках соответствует 00:57:00 ДМВ 18.09.2007.

Рис. 1.17. Сглаживание с помощью фильтра Калмана. Шаг по времени 200с. Момент = 0 на графиках соответствует 00:57:00 ДМВ 18.09.2007.

Рис. 1.18. Сглаживание с помощью фильтра Калмана. Шаг по времени 200с. Момент = 0 на графиках соответствует 00:57:00 ДМВ 18.09.2007.

Рис. 1.19. Сглаживание с помощью фильтра Калмана. Шаг по времени 200с. Момент = 0 на графиках соответствует 00:57:00 ДМВ 18.09.2007.

Рис. 1.20. Сглаживание с помощью фильтра Калмана. Шаг по времени 200с. Момент = 0 на графиках соответствует 00:57:00 ДМВ 18.09.2007.

Рис. 1.21. Сглаживание с помощью фильтра Калмана. Шаг по времени 200с. Момент = 0 на графиках соответствует 00:57:00 ДМВ 18.09.2007.

Определение вращательного движения спутника по данным измерений МПЗ и вектора угловой 2.1. Введение В главе 1 применялась методика, основанная на полных (динамических и кинематических) уравнениях вращательного движения спутника. Эти уравне­ ния содержат явные выражения моментов приложенных к спутнику внешних сил, и предназначены для описания неуправляемого движения. Такие выраже­ ния могут оказаться не достаточно точными, поэтому желательно проверить правильность реконструкции движения другими средствами. Проверку можно выполнить, обработав совместно измерения датчиков угловой скорости и борто­ вых магнитометров с помощью кинематической модели вращательного движе­ ния спутника. В основе этой модели лежат кинематические уравнения враща­ тельного движения для компонент кватерниона Q (см. п. 1.2). Входящие в эти уравнения компоненты угловой скорости строятся по данным измерений. Прин­ ципиальная возможность построения методики для определения вращательно­ го движения на основе только кинематических уравнений известна давно. В частности, в [36] с помощью математического моделирования показано, что из­ мерения напряженности МПЗ и угловой скорости КА позволяют определять движение последнего в реальном времени. В данной главе описаны два подхо­ да к реконструкции вращательного движения спутника: интегральная стати­ стическая методика, использующая метод наименьших квадратов, и методика, использующая фильтрацию Калмана. С помощью методики, основанной на ме­ тоде наименьших квадратов, было реконструировано вращательное движение тационной работе приведен обновленный вариант этой методики. В ее рамках данные измерений обоих типов, собранные на некотором отрезке времени, об­ рабатываются совместно. Эти данные сглаживаются тригонометрическими по­ линомами, которые подставляются в кинематические уравнения для компонент кватерниона, задающего ориентацию связанной со спутником системы коорди­ нат относительно гринвичской системы. Полученные таким образом уравнения представляют собой кинематическую модель вращательного движения спутни­ ка. Измерения МПЗ выбираются внутри отрезка времени, на котором определе­ ны эти уравнения. Решение кинематических уравнений вращательного движе­ ния, реконструирующее фактическое движение спутника, находится из условия наилучшего, в смысле метода наименьших квадратов, согласования данных из­ мерений вектора напряженности МПЗ с его расчетными значениями. Методика, использующая фильтрацию Калмана, разработана на перспективу. В дальней­ шем полет спутников научного назначения, создаваемых ФГУП "ЦСКБ-Про­ гресс", будет ориентированным (ориентация спутников солнечными батареями на Солнце будет поддерживаться двигателями маховиками или гиродинами), и методика мониторинга, основанная на кинематических уравнениях, станет основной. Для нее возникнет ситуация, упомянутая при описании задачи в гла­ ве 1.

Результаты второй главы опубликованы в работах [37–43].

2.2. Измерения угловой скорости на спутнике Фотон М- На последних трех спутниках типа были установлены трехос­ ные датчики угловой скорости системы управления движением. В случае КА — аппаратура DIMAC, которая предназначалась для измерения мик­ Фотон-М роускорений на борту спутника. Основными ее датчиками были два акселеро­ метра — низкочастотный и высокочастотный. Этой системой было проведено сеансов измерений, которые охватывали отрезки времени длиной 84 мин. Прак­ тически сразу после сеанса, данные измерений по телеметрическому каналу передавались на Землю.

Данные измерений угловой скорости интерпретировались в жестко связан­ Точка — центр масс ИСЗ, ось 1 параллельна продольной оси ИСЗ и на­ правлена от спускаемого аппарата к приборному отсеку.

Данные измерений, полученные во время одного сеанса, представляют со­ бой совокупность чисел где ( = 1, 2, 3) — значения компонент угловой скорости спутника в при­ борной системе координат в момент времени, = 12 с.

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

Аппроксимация данных (2.1) выполнялась с помощью дискретных рядов Фу­ рье [38–40] независимо для каждой векторной компоненты. Последовательности где, — коэффициенты и число одинаково для всех = 1, 2, 3. Это число должно не превосходить 1 и быть таким, чтобы выражения (2.2) позволяли с высокой точностью аппроксимировать на отрезке компоненты уг­ ловой скорости в уравнениях вращательного движения КА. Коэффициенты, находились методом наименьших квадратов. Полученные таким образом выра­ жения (2.2) иногда испытывают заметные высокочастотные колебания. Чтобы избавиться от них коэффициенты при старших гармониках корректировались с помощью специальных множителей [7, 8, 38–40] для ослабления присутствия в (2.2) высоких частот:

Здесь 1 — целая часть числа /2. Точность аппроксимации данных измере­ ний выражениями (2.2) характеризовалась стандартными отклонениями [45] В расчетах принималось 20 30, а среднеквадратичные ошибки равенств 2.3. Кинематическая модель движения спутника Спутник считаем твердым телом, геоцентрическое движение центра масс которого — кеплерово эллиптическое [46]. Элементы этого движения определя­ ются по данным траекторных измерений (ср. [7]). Для записи уравнений движе­ ния спутника относительно центра масс и соотношений, необходимых при обра­ ботке данных измерений, будем использовать введенную выше приборную си­ няя обозначения, введенные в главе 1, положение системы координат 1 относительно гринвичской системы координат будем задавать единичным ква­ тернионом Q.

Кинематические уравнения вращательного движения запишем в кватерни­ онной форме (ср. 1.2) Здесь — компоненты вектора абсолютной угловой скорости спутника в систе­ ме 1 2 3, () — выражения (2.2), аппроксимирующие измерения угловой скорости, — постоянные смещения в этих измерениях, — смещение шкалы времени аппаратуры DIMAC относительно шкалы времени системы управления движением спутника, и имеют такой же смысл, как и в (1.1). Параметры и считались неизвестными и определялись из обработки данных измерений МПЗ наряду с начальными условиями движения спутника.

2.4. Методика реконструкции вращательного движения спутника, по данным измерений его угловой скорости и вектора напряженности МПЗ Положим ( ) = +, +. Уравнения (2.3) рассматриваются при значениях и, удовлетворяющих условиям ( ),, ( ).

Следуя методу наименьших квадратов, аппроксимацией фактического дви­ жения спутника на отрезке ( ) будем считать решение системы (2.3), достав­ ляющее минимум функционалу [38–40] (ср. (1.4)) Здесь — постоянные смещения в псевдоизмерениях МПЗ, () вычисля­ лись по формуле (1.4), — число элементов множества ( ).

Отыскание решения системы (2.3), реконструирующего фактическое вра­ щательное движение спутника, состоит в минимизации функционала (2.4) по начальным условиям этого решения ( + ), и смещениям, при огра­ ничениях, наложенных условием нормировки (ср. (1.5)) Введем вектор R7 составленный из ( + ), и будем рассматривать функционал (2.4) как функцию (, ). Минимизация по и сводилась к вычислению функции в последовательности точек, сходящейся к пределу * = arg min 1 ( ) [38– 40]. Минимизация по при фиксированном выполнялась в два этапа (ср.

глава 1). На начальном этапе применялся метод Левенберга–Марквардта, далее метод Гаусса–Ньютона (ср. глава 1). На каждом шаге этих методов поправки к имеющимся значеним начальных условий решения системы (2.3) ищутся в виде (1.6). Расчет матрицы и свободного члена системы нормальных уравнений, возникающей на каждой итерации этих методов, осуществлялся по формулам, которые получены модификацией (1.3) Теперь 1, 2,..., 6 — обозначения величин 1, 2, 3, 1, 2, 3. Значения псевдопроизводных / определяются в процессе совместного интегриро­ вания уравнений и уравнений (2.3). Ненулевые значения производных / и начальных условий для / имеют вид Интегрирование этих уравнений выполнялось методом, описанным в разде­ ле 1.3.

Для нахождения начального приближения функционал (2.4) минимизи­ ровался на подмножестве допустимого множества значений параметров (, ), точки которого удовлетворяют условиям = 0, а ( + ) выражаются через параметр по формулам (1.3). На этом подмножестве функционал рассматри­ вался, как функция четырех параметров: и ( = 1, 2, 3). В таком случае начальные условия для уравнений в вариациях (2.7) имеют вид (1.3).

Выделение переменной из полного набора аргументов вызвано тем об­ стоятельством, что в процессе изменения множество ( ) может меняться.

В точках, где происходит это изменение функция, 1 ( ) недифференцируема.

По этой причине вычисление точек выполнялось без использования произ­ водной /.

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

Для оценки точности аппроксимации движения спутника применялся та­ кой же подход как и в разделе 1.3. Эта точность характеризовалась стандартны­ ми отклонениями ошибок аппроксимации псевдоизмерений МПЗ и найденных оценок и, полученными в рамках метода наименьших квадратов. Принятый способ расчета указанных характеристик соответствует следующим допущени­ ям [47]: 1) ошибки в псевдоизмерениях некоррелированы и имеют одинаковые дисперсии, 2) средние значения ошибок, относящихся к одной и той же компо­ ненте напряженности МПЗ, равны. Этот способ выбран из соображений удоб­ ства и вида функционала (2.4). Теоретико-вероятностные условия его адекват­ ности требуют дополнительного обсуждения.

Приведем расчетные формулы. Пусть (*, * ) = argmin (, ), — матри­ ца вычисленная в точке (*, * ). Тогда стандартное отклонение ошибок аппроксимации измерений находится по формуле ковариационная матрица оценок = *, = * равна 2 1. Стандарт­ ные отклонения этих оценок суть квадратные корни из соответствующих диагональных элементов ковариационной матрицы. Рассчитанные описанным способом стандартные отклонения оценок величин,, будем обозна­ 2.5. Примеры реконструкции неуправляемого движения методики из раздела 2.4 приведены в табл. 2.1 — 2.3 и на рис. 2.1 — 2.6. В табл. 2.1 для каждого имеющегося отрезка данных измерений угловой скоро­ сти приведены номер витка орбиты — В, день сентября 2007 г. — Д, скоррек­ тированные начальная — + и конечная точки — +, среднее квадра­ тическое отклонение аппроксимации измерений МПЗ и число включенных в обработку измерений на заключительном этапе минимизации функционала (2.4). В табл. 2.2 приведены значения углов,, и в момент времени +, смещений и. В табл. 2.3 приведены значения средних квадратических отклонений компонент вектора бесконечно малого поворота в момент времени +, смещения и. Единицы измерения величин, приведенных в табли­ Рисунки иллюстрируют найденное движение спутника и точность аппрок­ симаций измерений (1.3) функциями (). Время на рисунках измеряется в минутах от момента + до +. Все рисунки естественным образом раз­ биваются на две части — левую и правую. В левых частях рис. 2.1, 2.3 и 2. приведены графики зависимости от времени углов, и, а также график разности () = () 0 1 ( ), где 0 + 1 ( ) — линейная аппрокси­ мация функции (), построенная методом наименьших квадратов. В правых частях этих рисунков помещены графики выражений (), аппроксимирую­ щих данные измерений угловой скорости, и ломаные, проходящие через точки В левых частях рис. 2.2, 2.4 и 2.6 приведены графики функций () + (см. (2.4)) и ломаные, проходящие через точки (, ), ( ). Каждая ло­ маная и график аппроксимирующей ее функции () + изображены в еди­ ной системе координат. В выбранном масштабе ломаные и аппроксимирующие их графики совпадают, поэтому в правых частях рис. 2.2, 2.4 и 2.6 приведены Точность реконструкции реального движения спутника характеризуется стандартными отклонениями уточняемых параметров математической модели (табл. 2.3). Несмотря на то что, как уже отмечалось, обычные теоретико-веро­ ятностные допущения метода наименьших квадратов в данном случае не вы­ полнены, стандартные отклонения оказываются полезными характеристиками.

В частности, они позволяют сравнить точность определения различных пара­ метров. Анализ стандартных отклонений вместе анализом собственных векто­ ров матрицы (, ) (см. раздел 2.4), отвечающих ее нескольким минималь­ ным собственным числам, позволяет понять характер ошибок, возникающих при определении движения спутника. Анализ основан на формуле где и — собственные числа и ортонормированные собственные векторы матрицы. Из этой формулы следует, что наибольшие стандартные отклоне­ ния имеют те определяемые параметры, которым отвечают наибольшие компо­ ненты векторов /.

В качестве примера исследуем точность определение движения на данных, относящихся к витку 128 (2.5 и 2.6). Из приведенных в табл. 2.3 стандартных отклонений видно, что в указанных выше единицах интегрирования уравнений (2.3) наиболее точно определяется параметр, наименее точно определяются параметры 2, 3 : = 0.0001, 2 = 0.027, 3 = 0.032. Расположенные в порядке возрастания собственные числа матрицы в данном случае составля­ ют 1 = 10.2, 2 = 19.2, 3 = 64.9, 4 = 252.6, 5 = 359.7, 6 = 949.1, 7 = 14529. 1 = [ 0.0546, 0.0994, 0.0298, 0, 0496, 0.1029, 0.2675, 0.0019 ], 2 = [ 0.0458, 0.0111, 0.0326, 0.0061, 0.2096, 0.0697, 0.0000 ], 3 = [ 0.0601, 0.0679, 0.0267, 0.0636, 0.0065, 0.0488, 0.0020 ], 4 = [ 0.0495, 0.0350, 0.0129, 0.0046, 0.0067, 0.0002, 0.0074 ], 5 = [ 0.0003, 0.0172, 0.0483, 0.0034, 0.0092, 0.0077, 0.0011 ], 6 = [ 0.0086, 0.0136, 0.0067, 0.0272, 0.0022, 0.0019, 0.0006 ], 7 = [ 0.0007, 0.0006, 0.0004, 0.0000, 0.0000, 0.0000, 0.0082 ].

В этих расчетах (как и в программах минимизации функций и 1 ) единица измерения напряженности магнитного поля равнялась 50000. Компоненты векторов здесь упорядочены также, как в векторе (, ). Как видно из приве­ денных формул, пятая и шестая компонента вектора 1 / 1 и пятая компонента вектора 2 / 2, отвечающие параметрам 2 и 3, дают наибольший вклад в диагональные элементы матрицы 1. Минимальный вклад в эти элементы дают последние компоненты векторов /, отвечающие параметру.

Рис. 2.1. Движение Фотона реконструированное методом наименьших квадратов. Момент = 0 на графиках соответствует 14:11: ДМВ 15.09.2007.

Рис. 2.2. Аппроксимация магнитных измерений, полученная методом наименьших квадратов, = 1117, = 479. Момент = 0 на графиках соответствует 14:11:09 ДМВ 15.09.2007.

Рис. 2.3. Движение Фотона реконструированное методом наименьших квадратов. Момент = 0 на графиках соответствует 14:05: ДМВ 17.09.2007.

Рис. 2.4. Аппроксимация магнитных измерений, полученная методом наименьших квадратов, = 1092, = 544. Момент = 0 на графиках соответствует 14:05:48 ДМВ 17.09.2007.

Рис. 2.5. Движение Фотона реконструированное методом наименьших квадратов. Момент = 0 на графиках соответствует 12:27: ДМВ 22.09.2007.

Рис. 2.6. Аппроксимация магнитных измерений, полученная методом наименьших квадратов, = 1109, = 567. Момент = 0 на графиках соответствует 12:27:24 ДМВ 22.09.2007.

Таблица 2.1. Временные отрезки обработки с помощью кинематической модели Таблица 2.2. Оценки вектора параметров кинематической модели 1.5139 2.9213 1.1958 0.0509 0.0062 0. Таблица 2.3. Средние квадратические отклонения оценок параметров кинематической моде­ 17 0.0006 0.0012 0.0009 0.0012 0.0006 0.0010 0. 33 0.0009 0.0006 0.0010 0.0013 0.0009 0.0005 1. 49 0.0009 0.0008 0.0010 0.0011 0.0015 0.0017 0. 65 0.0009 0.0012 0.0006 0.0011 0.0014 0.0012 0. 81 0.0009 0.0013 0.0006 0.0010 0.0011 0.0011 0. 96 0.0009 0.0009 0.0007 0.0006 0.0011 0.0014 0. 113 0.0010 0.0011 0.0006 0.0009 0.0027 0.0024 0. 128 0.0012 0.0014 0.0008 0.0010 0.0027 0.0032 0. 144 0.0008 0.0010 0.0008 0.0007 0.0032 0.0029 0. 2.6. Фильтр Калмана Поскольку в дальнейшем полет спутников научного назначения будет ори­ ентированным (ориентация спутников солнечными батареями на Солнце будет поддерживаться двигателями маховиками или гиродинами), методика монито­ ринга, основанная на кинематических уравнениях, станет основной. Для нее возникнет ситуация, упомянутая при описании первой задачи. В связи с этим рассмотрим применение фильтра Калмана для кинематической модели.

При реализации фильтра, основанного на решении системы (1.13) компо­ ненты вектора состояния R10 механической системы представляют собой величины ( = 0, 1, 2, 3),, ( = 1, 2, 3). В отличие от пункта 2.4, смещения в измерениях компонент векторов напряженности МПЗ и угловой скорости не являются постоянными величинами. Теперь, эти смещения пред­ ставляют собой случайные величины.

Значения суть значения переменных вектора состояния в узлах времен­ являются серединами отрезков [, ] Функция () из (1.5) строится посредством интегрирования уравне­ ния (2.3), от точки 1 к точке. Проинтегрируем левую и правую часть уравнения (2.3) по времени на отрезке 1,, считая его решение известным.

Получим С учетом равномерности сетки { }, для гладких решений уравнения (2.3) име­ Последнее равенство представляет собой стандартную формулу трапеций [48], поэтому с погрешностью можно записать Из (2.10) следует Непосредственной проверкой, получим = 1, следовательно разностная схема (2.11) сохраняет норму кватерниона. Если () почти постоянна на отрез­ ке 1,, то формулу (2.11) можно использовать для вычисления значения кватерниона Q() в точке 1,, заменив на 1.

Применяя схему (2.11) для интегрирования системы уравнений вращатель­ ного движения (2.3), функцию (1 ) из (1.12) запишем в виде на полуинтервале [1, ) ошибки в измерениях компонент векторов МПЗ и угловой скорости соответственно. Следуя подходу принятому в пункте 2.4, сме­ щения и постоянны с точностью до случайных ошибок математической модели.

Формирование векторов измерений и вычисление функции () выпол­ няется следующим образом. Измерения напряженности магнитного поля (ср.

(1.3)) [, +1 ) — объединим в вектор. Для вычисления функции расчетных аналогов измерений, необходимо иметь значения кватерниона ори­ ентации Q( ) вычисленные в точке. Эти значения вычислялись по рекур­ рентной формуле (2.11) в предположении, что на полуинтервале [, +1 ) угло­ вая скорость была постоянна и равнялась ( ).

Матрицы Якоби функций () и () не имеют полного ранга, посколь­ ку компоненты кватерниона зависимы. Чтобы избежать вырождения, разность между компонентами фактического и прогнозируемого значений кватерниона ( ) записывалась в терминах вектора бесконечно малого поворота (см. 1.6).

Таким образом, функция (1.5) с высокой точностью можно записать в виде Здесь R9 — вектор составленный из компонент вектора бесконечно малого поворота и поправок, к прогнозируемым значениям и, — значение вектора поправок в узле. При использовании данной замены, вместо ковариаций компонент кватернионов используются ковариации компо­ нент векторов бесконечно малых поворотов. Прогноз матрицы ковариаций вычислялся по формуле (1.5), где теперь матрица Якоби = (^1 )/.

Матрицы,, и т.д. имеют порядок 15.

При минимизации функции (2.12), на каждом шаге решается система нор­ мальных уравнений Здесь — имеющееся приближенное значение оценки величины, которо­ му после решения системы уравнений в вариациях присваивается значение + ; R10 — поправка прогноза вектора состояния, отвечающая.

Когда норма поправки станет приемлемо малой, принимается = *.

Матрицы Якоби и вычислялись следующим образом. Матрица,=1 содержит следующие ненулевые элементы 3. Равенство нулю остальных элементов следует из предположений модели о независимости постоянных сдвигов в измерениях МПЗ и угловой скорости, а также теоремы о существовании и единственности решения задачи Коши.

Элементы ( )/ матрицы вычисляются интегрированием уравне­ ний в вариациях (1.8) с ненулевыми начальными условиями Для интегрирования системы дифференциальных уравнений вариациях при­ менялась разностная схема, аналогичная (2.11). Чтобы построить эту схему, запишем (1.8) в векторном виде где / = e (e — орт оси ) и / = 0 при всех остальных значениях. Интегрируя уравнение (2.14) на отрезке [1, ], получим Из (2.15) получим уравнение для вычисления ( )/ с погрешностью ( ) Векторное уравнение (2.16) удобно рассматривать, как систему линейных ал­ гебраических уравнений. Эта система решалась методом Гаусса с выбором ве­ дущего элемента по столбцу.

Как и в случае разностной схемы (2.11), при слабо меняющейся на от­ резке [1, ] угловой скорости (), формулу (2.16) можно использовать и для вычисления значения псевдопроизводных ()/ в произвольной точке Матрица =,=1 из (1.13) содержит следующие ненулевые эле­ менты Равенство нулю остальных элементов следует из предположений модели о неза­ висимости компонент постоянных сдвигов в измерениях МПЗ, а также теоремы о существовании и единственности решения задачи Коши.

Элеметы матрицы вычислялись посредством интегрирования уравне­ ний в вариациях (2.14) на отрезке, +1 с начальными условиями (2.6) в точке. Эти уравнения интегрировались при помощи модификации разност­ ной схемы (2.16), в предположении о постоянстве угловой скорости на отрезке интегрирования.

Точность аппроксимации измерений вектора МПЗ на каждом полуинтер­ вале [, +1 ) характеризуется величинами где — количество элементов множества ().

Оценка, матрица 01 и смещение шкалы времени магнитометров отно­ сительно шкалы времени датчиков угловой скорости находились методом наи­ меньших квадратов, посредством минимизации функционала (2.4). В обозначе­ ниях пункта 2.4 = *, 0 = 2.

2.7. Примеры реконструкции неуправляемого движения с помощью фильтра Калмана Продемонстрируем результаты реконструкций вращательного движения на трех отрезках времени, отвечающих 17, 49 и 128 виткам (см. табл. 2.1). Ри­ сунки организованы нижеописанным образом. Левые части рис. 2.7, 2.11, 2. содержат графики зависимости от времени углов, и, а также график разно­ сти () = () 0 1 ( ), где 0 + 1 ( ) — линейная аппроксимация функции (), построенная методом наименьших квадратов (ср. рис. 2.1, 2. и 2.5). В правых частях — ломаные, проходящие через точки (, + ( )).

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

В левых частях рис. 2.8, 2.12, 2.16 приведены графики ломаных, прохо­ дящие через точки (, ) и (, ( ) + ( )). В выбранном масштабе ломаные совпадают, поэтому в правых частях этих рисунков приведены лома­ На рис. 2.9, 2.13, 2.17 изображены графики величин, и. Графи­ ки демонстрируют сходимость фильтра Калмана. Эти величины — кусочно постоянные функции, на полуинтервалах [, +1 ) они сохраняют свои значе­ ния. На рис. 2.10, 2.14, 2.18 приведены графики стандартных отклонений оценок. Каждая компонента оценки представлена соответствующим стандартным отклонением, за исключением компонент кватерниона. Последние представле­ ны стандартными отклонениями компонент вектора бесконечно малого по­ ворота ( = 1, 2, 3). Это также кусочно постоянные функции с интервалами постоянства [, +1 ).

Рис. 2.7. Движение реконструированное c помощью фильтра Калмана. Момент = 0 на графиках соответствует 14:11: ДМВ 15.09.2007.

Рис. 2.8. Движение реконструированное c помощью фильтра Калмана. Момент = 0 на графиках соответствует 14:11: ДМВ 15.09.2007.

Рис. 2.9. Движение реконструированное c помощью фильтра Калмана. Момент = 0 на графиках соответствует 14:11: ДМВ 15.09.2007.

Рис. 2.10. Движение реконструированное c помощью фильтра Калмана. Момент = 0 на графиках соответствует 14:11: ДМВ 15.09.2007.

Рис. 2.11. Движение реконструированное c помощью фильтра Калмана. Момент = 0 на графиках соответствует 14:05: ДМВ 17.09.2007.

Рис. 2.12. Движение реконструированное c помощью фильтра Калмана. Момент = 0 на графиках соответствует 14:05: ДМВ 17.09.2007.

Рис. 2.13. Движение реконструированное c помощью фильтра Калмана. Момент = 0 на графиках соответствует 14:05: ДМВ 17.09.2007.

Рис. 2.14. Движение реконструированное c помощью фильтра Калмана. Момент = 0 на графиках соответствует 14:05: ДМВ 17.09.2007.

Рис. 2.15. Движение реконструированное c помощью фильтра Калмана. Момент = 0 на графиках соответствует 12:27: ДМВ 22.09.2007.

Рис. 2.16. Движение реконструированное c помощью фильтра Калмана. Момент = 0 на графиках соответствует 12:27: ДМВ 22.09.2007.

Рис. 2.17. Движение реконструированное c помощью фильтра Калмана. Момент = 0 на графиках соответствует 12:27: ДМВ 22.09.2007.

Рис. 2.18. Движение реконструированное c помощью фильтра Калмана. Момент = 0 на графиках соответствует 12:27: ДМВ 22.09.2007.

Проверка согласованности данных измерений магнитометров, установленных на борту ИСЗ 3.1. Введение Обработка данных измерений магнитного поля Земли, полученных на бор­ ту ИСЗ, обычно выполняется с использованием достаточно сложных матема­ тических моделей. Желательно провести предварительную проверку имеющих­ ся данных простыми средствами. Если измерения проводились одновременно несколькими магнитометрами, то в качестве такой проверки можно использо­ вать проверку геометрической согласованности их показаний. Если проверка оказывается успешной, то в результате удается оценить постоянные смещения в измерениях и матрицы перехода между собственными системами координат магнитометров. В данной главе описана методика проверки согласованности измерений двух магнитометров. Приведены примеры ее применения при обра­ ботке данных, полученных аппаратурой DIMAC (см. п. 1.3). Подобная провер­ ка данных измерений МПЗ, полученных аппаратурой "Мираж", на спутнике была приведена в статьях [49, 50]. Ниже описан один из возмож­ Фотон М- ных вариантов такой проверки, который использовался при обработке данных, полученных аппаратурой DIMAC (см. п. 1.3).

Аппаратура DIMAC имела четыре трехкомпонентных магнитометра, уста­ новленных в разных частях спускаемого аппарата и занумерованных числами 0, 1, 2 и 3. Измерения проводились в течение всего полета. Оцифровка пока­ заний всех магнитометров выполнялась для одних и те же моментов времени.

Интервалы между соседними измерениями не были постоянными. Их длина варьируется в пределах 1 12 с, а в среднем составляет около 5 с.

Результаты третьей главы опубликованы в работах [49, 50].

3.2. Методика проверки согласованности данных измерений бортовых магнитометров Пусть на борту искусственного спутника Земли установлены два магнито­ метра. Обозначим их I и II. Магнитометры расположены в разных местах, но измеряют одно и то же поле — МПЗ. Требуется проверить согласованность по­ лученных данных. При этом допускается, что измерения обоих магнитометров могут содержать постоянные смещения.

Точная постановка задачи состоит в следующем. На некоторой временной сетке 0 < 1 < · · · < заданы два набора компонент вектора напряженности магнитного поля H. Компоненты относятся к собственной системе координат магнитометра I, компоненты — к собственной системе координат магнито­ собой данные измерений магнитометров I и II соответственно, выполненные с малыми ошибками в момент времени. Если пренебречь ошибками, то при каждом величины и являются компонентами одного и того же векто­ ра H( ) и поэтому связаны определенными соотношениями. Эти соотношения, игнорируя случайные ошибки в данных измерений и учитывая ошибки систе­ матические (постоянные смещения), можно записать в виде Здесь — пересчитанные в систему координат магнитометра I постоянные сме­ щения в данных измерений, — элементы матрицы перехода =,=1 от системы координат магнитометра II к системе координат магнитометра I. При наличии случайных ошибок в данных измерений соотношения (3.1) становятся приближенными.

Если на отрезке 0 спутник совершает сложное вращательное движение, то в общем случае уравнения (3.1) достаточны для определения мат­ рицы и смещений уже при 5. Однако в реальных ситуациях и 0 велики. Учитывая приближенность соотношений (3.1) и большое зна­ чение, для отыскания матрицы и смещений удобно воспользоваться методом наименьших квадратов. Применение этого метода означает принятие следующей гипотезы: ошибки в соотношениях (3.1) некоррелированы, имеют нулевые математические ожидания и одинаковые дисперсии. Следуя методу наименьших квадратов, будем искать минимум выражения по величинам и при условии, что матрица ортогональна и имеет положительный определитель. Решение несколько более простой задачи, ко­ гда = 0 ( = 1, 2, 3) и требуется минимизировать выражение (3.2) только по элементам матрицы, хорошо известно (см., например, [51]). Незначительная модификация этого решения, позволяет выполнить полную минимизацию.

Выражение (3.2) с учетом ортогональности представим в виде Поставленная задача минимизации является задачей на условный экстремум — при минимизации необходимо учитывать условия ортогональности матри­ Здесь — символ Кронекера. Для решения задачи воспользуемся методом неопределенных множителей Лагранжа. Составим функцию где = — неопределенные множители Лагранжа и в выражении для опущены слагаемые, не зависящие от и. Условия безусловного минимума по величинам и имеют вид Из первой группы этих условий находим Подставив результат во вторую группу условий, получим Последние соотношения запишем в матричном виде Здесь — ортогональная, а — симметричная матрицы. Уравнение (3.4) ре­ шается следующим образом [51]. Рассмотрим сингулярное разложение матри­ цы [52] : =. Здесь и — ортогональные матрицы порядка 3, = diag (1, 2, 3 ), 1 2 3 0. Полагаем, что 3 > 0, т. е. матрица не вырождена. Введем матрицу = diag (1, 2, 3 ), = ±1 ( = 1, 2, 3), но выбор знаков пока не фиксируем. На основании сингулярного разложения Это наиболее общий вид матриц и с требуемыми свойствами.

Полученные формулы определяют несколько решений уравнения (3.4).

Выберем из них то, которое доставляет минимум и удовлетворяет условию det = 1. На решениях уравнения (3.4) = 2tr + 0, где 0 не зависит от. Простые преобразования дают Поскольку det = 1 2 3 det det, следует положить 1 = 2 = 1, 3 = det det. Окончательное выражение для матрицы имеет вид После того, как матрица найдена, смещения вычисляются по форму­ лам (3.3). Для вычисления сингулярного разложения матрицы используется соответствующая подпрограмма из [53], переписанная на C#.

Найденное решение обозначим, =. Оценим его точность. С этой целью линеаризуем задачу минимизации выражения (3.2) в окрестности точки минимума. Малые ошибки в задании ориентации системы координат маг­ нитометра II по отношению к системе координат магнитометра I будем описы­ вать в терминах вектора ее бесконечно малого поворота = (1, 2, 3 ). Ком­ поненты этого вектора будем указывать в системе координат магнитометра I.

Элементы произвольной ортогональной матрицы =,, можно в линейном приближении по представить в виде Здесь — символ Леви-Чивиты.

Соотношения (3.1), учитывая наличие в них ошибок, представим в виде Здесь — ошибки. Положим Вычтем последние соотношения из соотношений (3.6) и в полученных равен­ ствах перегруппируем члены. Будем иметь Подставим сюда соотношения (3.5), при этом в левой части изменим порядок суммирования. В результате придем к равенствам ния будем рассматривать как линейную задачу метода наименьших квад­ ратов для определения величин и. Это и есть упоминавшая­ ся выше линеаризация задачи минимизации выражения (3.2). Введем вектор = (1, 2, 3, 1, 2, 3 ), определим подходящим образом и запишем задачу (3.7) в виде. Ее решение имеет вид = ( )1, ковариационная матрица оценки равна По смыслу линеаризации в данном случае = 0, (^ ) (^ ) — зна­ чение выражения (3.2) в точке минимума. Обозначим это значение min и по­ ложим 2 = min /3( 1) в выражении для, получим формулу оценки ковариационной матрицы =,=1 в нелинейной задаче. Определенная по-новому величина — оценка стандартного отклонения ошибок выполнения ния оцениваемых параметров; в данном случае при = 1, 2, 3 это — стандартные отклонения оценок смещений, при = 4, 5, 6 — среднеквадратичные зна­ чения углов 1, 2, 3. Указанные среднеквадратичные значения обозначим 1, 2, 3 и будем использовать как характеристики точности оценки. Знание величин, и позволяет обоснованно судить о согласованности показаний магнитометров I и II.

Наряду с описанным подходом рассмотрим классический подход. Матрицу параметризуем углами, и, которые введем с помощью следующего условия. При совмещении с помощью параллельного переноса начал систем координат обоих магнитометров система магнитометра I переводится в систему магнитометра II тремя последовательными поворотами вокруг своих осей: 1) на угол вокруг оси 2; 2) на угол вокруг оси 3, преобразованной первым поворотом; 3) на угол вокруг дважды преобразованной оси 1, совпадающей с осью 1 системы координат магнитометра II. Элементы матрицы выражаются через введённые углы по формулам Ниже для удобства записи формул используются обозначения: = 2, = 3, Отыскание углов,, и смещений сводится к минимизации по этим величинам выражения (3.2). Минимизация выполнялась методом Гаусса-Нью­ тона [28]. На каждой итерации этого метода поправки,, уточняющие имеющиеся оценки и, определяется системой Использованные здесь величины, ( = 1, 2, 3) играют роль компонент угло­ вой скорости, отвечающей параметру. Они служат для представления про­ изводных элементов матрицы по этому параметру:

(обычная угловая скорость служит для представления производных по време­ ни). Компоненты, относятся к системе координат магнитометра I.

Пусть по-прежнему min — значение выражения (3.2) в точке минимума.

Матрицу системы (3.8), вычисленную в этой точке, обозначим. Тогда кова­ риационная матрица параметров, имеет вид Введенные выше векторы бесконечно малого поворота и вариации углов связаны соотношениями Отсюда следует формула =, где = diag (3, ), 3 — единичная матрица порядка 3, матрица =,,=1 вычислена в точке минимума выра­ жения (3.2). Эта формула позволяет найти среднеквадратичные значения углов по стандартным отклонениям углов.

3.3. Примеры проверки согласованности показаний магнитометров Методики предыдущего раздела проиллюстрируем результатами проверки согласованности данных измерений магнитометров 1, 2 и 3 аппаратуры DIMAC на Фотоне Проверка выполнялась обоими описанными способами. Началь­ ным приближением для метода Гаусса-Ньютона в классическом способе служи­ ли и углы,,, рассчитанные по матрице. Результаты применения разных способов совпали. Некоторые полученные результаты представлены в табл. 3.1 – 3.4 и на рис. 3.1 — 3.3.

В табл. 3.1 указаны интервалы времени, на которых проводилась провер­ ка согласования данных. Первый столбец таблицы содержит номер интервала, второй и третий столбцы — декретное московское время (ДМВ) его начальной и конечной точек. Интервалы 1 и 12 относятся к началу и концу измерений соответственно и имеют длину около 10 часов. Измерения МПЗ, относящиеся к 18.09.2007 и 24.09.2007, из-за наличия в них лакуны разбиты на два интер­ вала. Каждый из остальных интервалов охватывает примерно сутки. Табл. 3. — 3.4 демонстрируют согласованность данных измерений различных пар маг­ нитометров. Табл. 3.2 составлена для пары (1, 2), табл. 3.3 — для пары (1, 3), табл. 3.4 — для пары (2, 3). В каждой перечисленной паре первый магнитометр выступал в роли магнитометра I, второй — в роли магнитометра II. Структура табл. 3.2 — 3.4 одинакова. В них для интервалов из табл. 3.1 приведены оценки параметров согласования,, и, стандартные отклонения этих парамет­ ров,, и, а также стандартное отклонение ошибок согласования данных. Размерности перечисленных величин: [ ] = [ ] = [] = = 105 Э, Таблица 3.1. Временные интервалы сравнения показаний магнитометров Таблица 3.2. Результаты сравнения показаний магнитометров 1 и 2, [] = [ ] = [ ] = нТ ( = 1, 2, 3) Таблица 3.3. Результаты сравнения показаний магнитометров 1 и 3, [] = [ ] = [ ] = нТ ( = 1, 2, 3) Таблица 3.4. Результаты сравнения показаний магнитометров 2 и 3, [] = [ ] = [ ] = нТ ( = 1, 2, 3) Рис. 3.1 — 3.3 иллюстрируют достигнутое согласие показаний магнитомет­ ров на некоторых интервалах из табл. 3.1. Момент времени = 0 на этих ри­ сунках соответствует моменту включения аппаратуры DIMAC — 14:10:20 ДМВ 14.09.2007. Левые части рисунков иллюстрируют согласие данных двух магни­ тометров на интервале 1, средняя — на интервале 7, правая — на интервале 13.

На рисунках в каждой системе координат изображены два графика. Графики практически сливаются. Один из них представляет собой ломаную с вершина­ ми в точках,, = 0, 1,...,, ординаты которых определены левыми частями формулы (3.1). Каждое звено ломаной соединяет две точки с индекса­ ми, отличающимися на 1. Другой график представляет собой аналогичную ломаную, ординаты вершин которой определены правыми частями той же фор­ мулы.

Рис. 3.4 иллюстрирует сравнение данных магнитометров 0 и 1. Его структу­ ра идентична рис. 3.1 — 3.3. Выделение рис. 3.4 обусловлено несколько худшей согласованностью измерений магнитометра 0 с измерениями остальных.

Рис. 3.1. Сравнение данных измерений магнитометров 1 и 2 на интервалах 1, 7 и Рис. 3.2. Сравнение данных измерений магнитометров 1 и 3 на интервалах 1, 7 и Рис. 3.3. Сравнение данных измерений магнитометров 2 и 3 на интервалах 1, 7 и Рис. 3.4. Сравнение данных измерений магнитометров 1 и 2 на интервалах 1, 7 и Судя по таблицам и рисункам, показания магнитометров 1, 2 и 3 хорошо согласованы между собой и являются достаточно точными измерениями МПЗ.

Большие значения смещений ( = 1, 2, 3) в табл. 3.2 — 3.4 можно объяснить наличием на борту спутника большого количества проводов с током. Показания магнитометра 0 несколько отличаются от показаний остальных магнитометров.

Возможно, различия в показаниях объясняются неисправностью магнитометра 0; возможно, показания магнитометра 0 искажены влиянием другого магнит­ ного поля. Магнитометры аппаратуры DIMAC, размещались в разных местах спускаемого аппарата, а на спутниках Фотон достаточно много приборов, гене­ рирующих локальные магнитные поля.

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

Учитывая это обстоятельство и принимая во внимание большие значения, при обработке магнитных измерений, выполненных магнитометрами 1, 2 и 3, приближенные равенства 0, 0, 0 были заменены соответствующи­ ми точными равенствами. Конкретные значения углов определялись выбором пары магнитометров. Измерениями компонент магнитного поля, полученными аппаратурой DIMAC в некоторый момент времени, считались величины где, и — показания магнитометров 1, 2 и 3 в их собственных систе­ мах координат, полученные в тот же момент. Компоненты относятся к систе­ ме координат магнитометра 1. Интерпретированные таким способом измерения магнитометров аппаратуры DIMAC, позволили выполнить реконструкцию вра­ Основные результаты работы:

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