Космические исследования, 2019, T. 57, № 1, стр. 61-73

1. МАЛЫЙ КОСМИЧЕСКИЙ АППАРАТ АИСТ-2Д

Малый космический аппарат (МКА) Аист-2Д был разработан АО РКЦ “Прогресс” и выведен на орбиту 28.IV.2016 г. в 05:01:21 декретного московского времени (ДМВ) с космодрома “Восточный” ракетой-носителем “Союз-2” этапа 1а с использованием блока выведения “Волга”. Блоком выведения была сформирована целевая орбита МКА с параметрами: драконический период 5653.55 с, наклонение 97.27°, средняя высота 491.55 км, эксцентриситет 0.0004, угол между плоскостью орбиты и направлением на среднее Солнце 348.75°.

На борту МКА установлены: целевая оптико-электронная аппаратура дистанционного зондирования Земли “Аврора” среднего разрешения и комплекс научной аппаратуры (НА), включающий аппаратуру КМУ-1 (Компенсатор Микро Ускорений). Аппаратура КМУ-1 является результатом совместной разработки Института космического приборостроения Самарского университета им. С.П. Королёва и РКЦ “Прогресс”.

Функционирование НА в орбитальном полете поддерживается бортовыми обеспечивающими системами: системой контроля и управления, системой энергопитания, системой обеспечения теплового режима, системой управления движением, бортовой кабельной сетью. Передача данных от НА осуществляется по высокоскоростной радиолинии на наземные средства управления, получения и обработки информации (НСУ ПОИ) РКЦ “Прогресс”.

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

Включение основного комплекта НА КМУ-1 было осуществлено в 09:40 ДМВ 9.VI.2016 г. Основной комплект был выключен в 13:30 12.VIII.2016 г. Резервный комплект был включен в 14:30 ДМВ 12.VIII.2016 г. и работает до настоящего времени. Поступающая от НА информация ежесуточно обрабатывается программными средствами автоматизированного рабочего места КМУ-1 для контроля оперативных параметров, характеризующих работоспособность НА и ее отдельных структурных элементов.

Состав данных, переданных НА, позволяет решать задачу реконструкции вращательного движения МКА с использованием методического и программного обеспечения, разработанного в ИПМ им. М.В. Келдыша РАН. Ранее подобное программное обеспечение применялось для реконструкции вращательного движения спутников Бион М-1 и Фотон М-4 [1–3]. Его применение в случае МКА Аист-2Д потребовало доработки для определения матриц перехода от собственных систем координат магнитометров к приборной системе координат спутника, в которой задаются измерения угловой скорости. Ниже приводятся первые результаты реконструкция вращательного движения МКА (для удобства называемого спутником), полученные по данным НА КМУ-1.

2. ТЕСТИРОВАНИЕ МАГНИТНЫХ ИЗМЕРЕНИЙ

Магнитометры, установленные на спутнике, будем называть магнитометр 1 и магнитометр 2. Магнитные измерения проводились на отрезках полета длиной около 10 ч. Оцифровка показаний магнитометров выполнялась в единые моменты времени с шагом 21–22 с. Компоненты измеряемых векторов напряженности магнитного поля выдавались в собственных системах координат магнитометров. Собственные системы координат магнитометров 1 и 2 обозначим соответственно ${{y}_{1}}{{y}_{2}}{{y}_{3}}$ и ${{z}_{1}}{{z}_{2}}{{z}_{3}}.$ Система ${{z}_{1}}{{z}_{2}}{{z}_{3}}$ – это штатная система координат магнитометра 2, система ${{y}_{1}}{{y}_{2}}{{y}_{3}}$ переводится в штатную систему координат магнитометра 1 преобразованием ${{y}_{1}} \leftrightarrow {{y}_{2}},$ ${{y}_{3}} \to - {{y}_{3}}.$ Оси системы ${{y}_{1}}{{y}_{2}}{{y}_{3}}$ составляют углы менее 5° с одноименными осями приборной системы координат спутника ${{x}_{1}}{{x}_{2}}{{x}_{3}}$ и острые углы с одноименными осями собственной системы координат магнитометра 2.

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

Здесь $g_{i}^{{(n)}}$ и $h_{i}^{{(n)}}$ – результаты измерений i-х компонент магнитного поля магнитометрами 1 и 2 в их собственных системах координат в момент времени ${{t}_{n}},$ ${{t}_{1}} < {{t}_{2}} < \ldots < {{t}_{N}}.$ Тестирование измерений магнитометров выполнялось двумя способами.

Первый способ основан на предположении, что магнитометры измеряли в основном МПЗ. С учетом большого числа электрических цепей на борту спутника это предположение выглядит несколько наивным, однако, опыт магнитных измерений на двух предыдущих Аистах и ряде других спутников [1–6] позволял надеяться, что оно окажется адекватным. Этот способ состоял в сравнении модуля напряженности измеренного поля с модулем напряженности МПЗ, рассчитанной по аналитической модели IGRF. Сравнение модулей не требует знания ориентации магнитометра, надо знать только орбитальное движение спутника. Такой подход позволяет проверить данные измерений каждого магнитометра, используя минимальную дополнительную информацию.

Рассмотрим магнитометр 1. Будем считать, что его измерения (1) содержат две систематические ошибки. Одна ошибка вызвана неточным заданием масштабирующего коэффициента при оцифровке, другая – наличием постоянных смещений в измерениях. Правильные показания магнитометра в момент ${{t}_{n}}$ должны были бы быть $\bar {g}_{i}^{{(n)}} = {{\kappa }_{1}}g_{i}^{{(n)}} - {{a}_{i}}$ $(i = 1,2,3).$ Значения величин ${{\kappa }_{1}}$ и ${{a}_{i}}$ находились методом наименьших квадратов из условия минимума выражения

где ${{H}^{{(n)}}}$ – расчетные значения модуля напряженности МПЗ на спутнике в момент времени ${{t}_{n}}.$ Величины ${{H}^{{(n)}}}$ рассчитываются с использованием аналитической модели IGRF, модели орбитального движения SGP4 [7] и двухстрочных элементов орбиты спутника. Минимизация $\Psi $ выполнялась методом Гаусса–Ньютона [8]. Найденное решение характеризовалось стандартным отклонением ${{\sigma }_{H}} = \sqrt {{{{{\Psi }_{{\min }}}} \mathord{\left/ {\vphantom {{{{\Psi }_{{\min }}}} {(N - 4)}}} \right. \kern-0em} {(N - 4)}}} $ ошибок согласования величин ${{H}^{{(n)}}}$ и а также стандартными отклонениями определяемых параметров – квадратными корнями из соответствующих диагональных элементов матрицы $\sigma _{H}^{2}{{F}^{{ - 1}}},$ где F – матрица системы нормальных уравнений, решаемых при реализации метода Гаусса–Ньютона. Матрица F вычисляется в точке минимума $\Psi $ и приблизительно равна матрице квадратичной формы ${{{{d}^{2}}\Psi } \mathord{\left/ {\vphantom {{{{d}^{2}}\Psi } 2}} \right. \kern-0em} 2}$ в этой точке. Величины $\bar {g}_{i}^{{(n)}} = {{\kappa }_{1}}g_{i}^{{(n)}} - {{a}_{i}}$ с найденными ${{\kappa }_{1}}$ и ${{a}_{i}}$ будем считать скорректированными измерениями магнитометра 1.

Аналогичным образом находились скорректированные измерения второго магнитометра $\bar {h}_{i}^{{(n)}} = {{\kappa }_{2}}h_{i}^{{(n)}} - {{b}_{i}}$ $(i = 1,2,3),$ где корректирующие параметры ${{\kappa }_{2}}$ и ${{b}_{i}}$ находились из условия минимума выражения

Найденное решение характеризовалось стандартным отклонением ${{\sigma }_{H}}$ ошибок согласования величин ${{H}_{n}}$ и

а также стандартными отклонениями параметров ${{\kappa }_{2}}$ и ${{b}_{i}}.$

Приведем результаты коррекции измерений на трех отрезках данных, перечисленных в табл. 1. Здесь указаны начальные точки ${{t}_{1}}$ этих отрезков, числа N точек с измерениями и длины отрезков ${{t}_{N}} - {{t}_{1}}.$ Корректирующие параметры ${{\kappa }_{1}},$ ${{\kappa }_{2}},$ ${{a}_{i}},$ ${{b}_{i}}$ для этих отрезков и оценки точности их определения приведены в табл. 2, 3. Стандартные отклонения параметров указаны в скобках рядом с их значениями. Результаты сравнения модулей скорректированных измерений магнитного поля и расчетных значений модуля МПЗ иллюстрируются графиками на рис. 1. Рис. 1а построен для магнитометра 1. Верхние графики на этом рисунке представляют собой ломаные. Здесь в каждой системе координат изображены две ломаные. Звенья одной из них соединяют соседние по времени точки $({{t}_{n}},{{g}^{{(n)}}}),$ звенья другой ломаной аналогичным образом соединяют точки $({{t}_{n}},{{H}^{{(n)}}}),$ $n = 1,2, \ldots ,N.$ Эти ломаные практически совпадают, поэтому для наглядности на нижних графиках изображены ломаные, проходящие через точки $({{t}_{n}},{{g}^{{(n)}}} - {{H}^{{(n)}}}).$ Рис. 1б аналогичным образом характеризует согласие величин ${{h}^{{(n)}}}$ и ${{H}^{{(n)}}}.$

Данные в табл. 2, 3 и графики на рис. 1 показывают, что измерения обоих магнитометров достаточно точно согласуются со значениями МПЗ из модели IGRF. При этом согласно значениям ${{\sigma }_{H}}$ в табл. 2, 3 магнитометр 1 оказался примерно в полтора раза точнее магнитометра 2. Однако исключение измерений магнитометра 2 из обработки было бы неосмотрительным. Как будет показано ниже, линейная комбинация измерений обоих магнитометров с правильно подобранными весами позволит точнее реконструировать движение спутника, чем использование показаний одного только магнитометра 1.

Из найденных корректирующих параметров ниже используются только параметры ${{\kappa }_{1}},$ ${{\kappa }_{2}}.$ Во-первых, они определены с малой относительной ошибкой. Во-вторых, на последующих этапах обработки магнитных измерений используются более адекватные этим этапам способы оценки смещений. После того как параметры ${{\kappa }_{1}},$ ${{\kappa }_{2}}$ найдены, величины $g_{i}^{{(n)}}$ и $h_{i}^{{(n)}}$ в (1) заменялись величинами $\kappa {}_{1}g_{i}^{{(n)}}$ и ${{\kappa }_{2}}h_{i}^{{(n)}}.$ Чтобы не вводить новых обозначений, последние величины будем обозначать $g_{i}^{{(n)}},$ $h_{i}^{{(n)}}$ и именно их считать измерениями магнитного поля.

Второй способ тестирования магнитных измерений используется в случае, когда оба магнитометра прошли первый тест. Это способ основан на следующих предположениях: 1) магнитометры измеряют одно и то же переменное поле, 2) постоянные составляющие поля в точках расположения магнитометров различны (или измерения магнитометров содержат постоянные смещения). В случае точной реализации этих предположений измерения магнитометра 2, пересчитанные в систему координат магнитометра 1, отличались бы от измерений последнего на постоянный вектор смещения. Из-за разного рода ошибок точного постоянства вектора смещения достичь невозможно, однако его вариации во времени должны быть малыми. Согласующие параметры – матрица перехода между системами координат магнитометров и постоянный вектор смещения – определялись методом наименьших квадратов. Методика определения этих параметров описана в [9] и состоит в следующем.

С высокой точностью должны выполняться соотношения

где ${{d}_{i}}$ – компоненты постоянного вектора смещения в системе координат ${{y}_{1}}{{y}_{2}}{{y}_{3}},$ ${{c}_{{ij}}}$ – элементы матрицы перехода $C$ от системы координат ${{z}_{1}}{{z}_{2}}{{z}_{3}}$ к системе координат ${{y}_{1}}{{y}_{2}}{{y}_{3}}$ (${{c}_{{ij}}}$ – косинус угла между осями ${{y}_{i}}$ и ${{z}_{j}}$). Если на отрезке ${{t}_{1}} \leqslant t \leqslant {{t}_{N}}$ вектор магнитного поля совершал сложное вращательное движение относительно спутника, то для отыскания матрицы C и смещений ${{d}_{i}}$ удобно воспользоваться методом наименьших квадратов. Применение этого метода означает принятие следующей гипотезы: ошибки в соотношениях (2) некоррелированы, имеют нулевые математические ожидания и одинаковые дисперсии.

Следуя методу наименьших квадратов, ищем минимум выражения

по величинам ${{c}_{{ij}}}$ и ${{d}_{i}}$ при условии, что матрица $C$ ортогональна и ее определитель равен 1. Эта задача решается с помощью стандартных процедур вычислительной линейной алгебры. Обозначим решение этой задачи $C^\circ ,$ $d^\circ = {{(d_{1}^{^\circ },d_{2}^{^\circ },d_{3}^{^\circ })}^{{\text{T}}}}.$ В линейном приближении окрестность решения параметризируется независимыми параметрами ${{\xi }_{i}},{{\theta }_{i}}$ $(i = 1,2,3)\,:$ ${{d}_{i}} = d_{i}^{^\circ } + {{\xi }_{i}},$ $C = {{E}_{\theta }}C^\circ ,$ где

Величины ${{\theta }_{i}}$ представляют собой компоненты вектора бесконечно малого поворота, характеризующего отличие матрицы C от ее оценки $C^\circ .$ Эти компоненты относятся к системе координат ${{y}_{1}}{{y}_{2}}{{y}_{3}}.$ Согласно методу наименьших квадратов параметры ${{\xi }_{i}},\,\,{{\theta }_{i}}$ образуют случайный вектор $\eta \in {{R}^{6}}$ с нулевым математическим ожиданием. Ковариационная матрица ${{K}_{\eta }}$ этого вектора выражается через матрицу P системы нормальных уравнений, получающейся линеаризацией исходной задачи по η в точке минимума выражения Z, и значение ${{Z}_{{\min }}}{\text{:}}$

Здесь ${{\sigma }_{0}}$ – стандартное отклонение ошибок в соотношениях (2). Стандартные отклонения величин ${{\xi }_{i}}$ и ${{\theta }_{i}}$ – квадратные корни из соответствующих диагональных элементов матрицы ${{K}_{\eta }}$ – обозначим ${{\sigma }_{{di}}}$ и ${{\sigma }_{{\theta i}}}.$

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

Оценки точности этого решения и аналогичных решений для интервалов 1 и 3 приведены в табл. 4.

Наглядное представление о результатах сопоставления дает рис. 2. На рис. 2а в каждой системе координат изображены две ломаные. Звенья одной из них последовательно соединяют точки $\left( {{{t}_{n}},g_{i}^{{(n)}}} \right),$ другая проходит через точки $\left( {{{t}_{n}},\hat {g}_{i}^{{(n)}}} \right),$ $n = 1,2, \ldots N.$ На вид эти ломаные очень близки. Чтобы наглядно показать их различие, на рис. 2б изображены разности этих ломаных, т.е. ломаные с вершинами в точках $\left( {{{t}_{n}},\Delta g_{i}^{{(n)}}} \right),$ где $\Delta g_{i}^{{(n)}} = g_{i}^{{(n)}} - \hat {g}_{i}^{{(n)}}.$ Рисунок показывает, что различие между показаниями магнитометров достаточно велико.

Если из данных (1) сформировать взвешенную напряженность поля в системе ${{y}_{1}}{{y}_{2}}{{y}_{3}}\,:$

то полученные значения будут заметно лучше согласовываться со значениями IGRF в рамках первого теста, чем скорректированные данные обоих магнитометров, взятые по отдельности. Подставим величины $\tilde {g}_{i}^{{(n)}}$ вместо $g_{i}^{{(n)}}$ в выражение, задающее функцию $\Psi $, и положим

Для интервала 1 из табл. 1 при $\rho = 2$ получим $\sigma _{H}^{'} = 671\gamma ,$ для интервала 2 при $\rho = 1.7$ получим $\sigma _{H}^{'} = 673\gamma ,$ в случае интервала 3 при $\rho = 1.8$ будем иметь $\sigma _{H}^{'} = 654\gamma .$ Эти значения заметно лучше значений ${{\sigma }_{H}}$ в табл. 2, 3.

Для решения задачи реконструкции движения спутника матрицу $C$ параметризируем углами γ, α и β, которые введем посредством следующего условия. Система координат ${{y}_{1}}{{y}_{2}}{{y}_{3}}$ может быть переведена в систему ${{z}_{1}}{{z}_{2}}{{z}_{3}}$ тремя последовательными поворотами: 1) на угол α вокруг оси ${{y}_{2}};$ 2) на угол β вокруг оси ${{y}_{3}},$ получившейся после первого поворота; 3) на угол γ вокруг оси ${{y}_{1}},$ получившейся после первых двух поворотов и совпадающей с осью ${{z}_{1}}.$ Элементы матрицы $C$ выражаются через введенные углы с помощью формул

Углы, параметризирующие указанные выше матрицы $C^\circ ,$ приведены в табл. 5. Здесь стандартные отклонения ${{\sigma }_{\gamma }},$ ${{\sigma }_{\alpha }}$ и ${{\sigma }_{\beta }}$ (указаны в скобках рядом с оценками углов γ, α и β) вычислены по стандартным отклонениям ${{\sigma }_{{\theta {\kern 1pt} i}}}$ на основании формул θ₁ = c₁₁dγ + $ + \sin \alpha d\beta ,$ ${{\theta }_{2}} = {{c}_{{21}}}d\gamma + d\alpha ,$ ${{\theta }_{3}} = {{c}_{{31}}}d\gamma + \cos \alpha d\beta .$

Графики, иллюстрирующие тестирование магнитометров на интервалах 1 и 3, приведены в [10].

3. МЕТОДИКА РЕКОНСТРУКЦИИ ВРАЩАТЕЛЬНОГО ДВИЖЕНИЯ СПУТНИКА

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

Для описания вращательного движения спутника используются две правые декартовы системы координат: приборная система ${{x}_{1}}{{x}_{2}}{{x}_{3}},$ в которой интерпретируются измерения угловой скорости, и система ${{Y}_{1}}{{Y}_{2}}{{Y}_{3}}$, в которой задаются двухстрочные элементы. Последняя переводится в гринвичскую систему координат поворотом по часовой стрелке вокруг оси ${{Y}_{3}}$ на угол, равный среднему звездному времени. Таким образом, система ${{Y}_{1}}{{Y}_{2}}{{Y}_{3}}$ весьма близка ко второй геоэкваториальной системе координат эпохи даты. В частности, плоскость ${{Y}_{1}}{{Y}_{2}}$ совпадает с плоскостью экватора, ось ${{Y}_{3}}$ – ось вращения Земли. Систему ${{Y}_{1}}{{Y}_{2}}{{Y}_{3}}$ считаем инерциальной.

Положение приборной системы ${{x}_{1}}{{x}_{2}}{{x}_{3}}$ относительно системы ${{Y}_{1}}{{Y}_{2}}{{Y}_{3}}$ зададим нормированным кватернионом $Q = ({{q}_{0}},{{q}_{1}},{{q}_{2}},{{q}_{3}}),$ $q_{0}^{2}\, + \,q_{1}^{2}\, + \,q_{2}^{2}\, + \,q_{3}^{2}\, = \,1.$ Кватернионная формула перехода между этими системами координат имеет вид

Матрицу перехода от приборной системы к системе ${{Y}_{1}}{{Y}_{2}}{{Y}_{3}}$ обозначим $A = \left\| {{{a}_{{ij}}}} \right\|_{{i,j = 1}}^{3},$ где ${{a}_{{ij}}}$ – косинус угла между осями ${{Y}_{i}}$ и ${{x}_{j}}.$ Элементы этой матрицы выражаются через компоненты Q с помощью формул

Орбитальное движение спутника на представляющем интерес временном интервале описывается с помощью модели SGP4. Параметры модели определяются по данным двухстрочных элементов. Математическая модель вращательного движения спутника построена на основе кинематических уравнений, которым удовлетворяет кватернион Q. Эти уравнения имеют вид

Здесь точка над символом означает дифференцирование по времени, ${{\omega }_{1}},$ ${{\omega }_{2}}$ и ${{\omega }_{3}}$ – компоненты абсолютной угловой скорости спутника в приборной системе координат.

Измерения величин ${{\omega }_{i}}$ выполнялись непрерывно в течение практически всего полета. Данные этих измерений представлены в цифровом виде на единой временной сетке с постоянным шагом 12 с. По этим данным построены непрерывные кусочно-линейные функции ${{\Omega }_{i}}(t)$ $(i = 1,\;2,\;3)$, которые в узлах сетки совпадают с данными измерений величин ${{\omega }_{i}},$ а в промежутках между узлами линейны. Примеры графиков функций ${{\Omega }_{i}}(t)$ приведены на рис. 3. Однако в уравнения (1) подставлялись не эти функции, а выражения

где ${{\chi }_{i}}$ – постоянные смещения в измерениях. Уравнения (4), (5) использовались на интервалах времени не более 12 ч, причем смещения ${{\chi }_{i}}$ служили уточняемыми параметрами. В такой ситуации предположение о постоянстве этих величин оказалось вполне адекватным.

Опишем применение уравнений (4), (5) для обработки магнитных измерений, выполненных внутри отрезка $I = \{ t\,:{{t}_{a}} \leqslant t \leqslant {{t}_{b}}\} ,$ на котором построены функции ${{\Omega }_{i}}(t).$ Узлы временной сетки с магнитными измерениями, как правило, не совпадают с узлами сетки, на которой заданы функции ${{\Omega }_{i}}(t).$ Пусть для определенности обрабатываются измерения магнитометра 1 из (1), причем $[{{t}_{1}},{{t}_{N}}] \subset ({{t}_{a}},{{t}_{b}}).$ Матрицу перехода от системы ${{x}_{1}}{{x}_{2}}{{x}_{3}}$ к системе ${{y}_{1}}{{y}_{2}}{{y}_{3}}$ обозначим $\left\| {{{b}_{{ij}}}} \right\|_{{i,j = 1}}^{3},$ где ${{b}_{{ij}}}$ – косинус угла между осями ${{y}_{i}}$ и ${{x}_{j}}.$ Эту матрицу параметризируем углами ${{\gamma }_{1}},$ ${{\alpha }_{1}}$ и ${{\beta }_{1}}$, по смыслу и формулам, принятым в п. 2 при параметризации матрицы C углами γ, α и $\beta $.

Следуя методу наименьших квадратов, реконструкцией фактического движения спутника на отрезке I будем считать решение уравнений (4), (5), доставляющее минимум функционалу

Здесь ${{\Delta }_{i}}$ – постоянные смещения в измерениях магнитометра 1, ${{H}_{i}}(t)$ – расчетные значения компонент напряженности МПЗ в системе координат ${{Y}_{1}}{{Y}_{2}}{{Y}_{3}}$ в момент времени t. Функции ${{H}_{i}}(t)$ строятся с использованием модели орбитального движения SGP4 и аналитической модели МПЗ IGRF. Функционал (6) получен из функционала метода наименьших квадратов, возникающего при уравнивании соотношений $g_{i}^{{(n)}}\, \approx \,{{\hat {h}}_{i}}({{t}_{n}})\, + \,{{\Delta }_{i}}$ $(i\, = \,1,2,3;$ $n\, = \,1,2, \ldots ,N),$ посредством предварительной минимизации этого последнего функционала по смещениям ${{\Delta }_{i}}.$ Минимизация $\Phi $ выполняется по начальным условиям решения ${{q}_{j}}({{t}_{a}})$ $(j = 0,1,2,3),$ параметрам ${{\chi }_{i}}$ $(i = 1,2,3)$ и углам ${{\gamma }_{1}},$ ${{\alpha }_{1}}$ и ${{\beta }_{1}}.$ При этом учитывается условие нормировки

Для простоты письма все уточняемые величины объединим в один вектор $z \in {{R}^{{10}}}.$ В принятых обозначениях $\Phi = \Phi (z),$ $z* = \arg \min \Phi (z)$ – искомая оценка вектора z. Минимизация $\Phi (z)$ выполнялась в несколько этапов разными методами. Начальные этапы, на которых ${{\gamma }_{1}} = {{\alpha }_{1}} = {{\beta }_{1}} = 0,$ описаны в [1, 2]. Здесь ограничимся описанием заключительного этапа, когда эти углы уточняются.

Описание минимизации функции Φ по z начнем с заключительного этапа, на котором применялся метод Гаусса–Ньютона [8]. На каждой итерации этого метода поправки $\Delta {{q}_{j}}({{t}_{a}})$ к имеющимся значениям ${{q}_{j}}({{t}_{a}})$ ищутся в виде (ср. уравнения (4))

Параметры ${{\theta }_{i}}$ суть компоненты вектора бесконечно малого поворота, задающего изменение ориентации спутника в окрестности положения $Q({{t}_{a}}).$ Эти параметры и поправки $\Delta {{\chi }_{i}}$ находятся из системы нормальных уравнений с матрицей $\left\| {{{G}_{{ij}}}} \right\|_{{i,j = 1}}^{9}$ и правой частью $\left\| {{{D}_{i}}} \right\|_{{i = 1}}^{9}\,:$

Здесь ${{p}_{1}},{{p}_{2}}, \ldots {{p}_{9}}$ – обозначения величин θ₁, θ₂, θ₃, ${{\chi }_{1}},$ ${{\chi }_{2}},{{\chi }_{3}},{{\gamma }_{1}},{{\alpha }_{1}},{{\beta }_{1}}$ в указанном порядке, $\partial {{\hat {\varphi }}_{l}}{{(t)} \mathord{\left/ {\vphantom {{(t)} {\partial {{p}_{i}}}}} \right. \kern-0em} {\partial {{p}_{i}}}}$ – псевдопроизводные, служащие для представления истинных производных. Псевдопроизводная – это не частная производная некоторой функции по какому-то параметру. Запись ее в виде частной производной используется лишь для удобства. Такую запись следует воспринимать как единый символ с двумя индексами. Псевдопроизводная – это вектор, являющийся аналогом угловой скорости. В кинематике твердого тела угловая скорость служит для расчета производных по времени, а псевдопроизводная используется для расчета производных по параметру. В обозначении $\partial {{\hat {\varphi }}_{l}}{{(t)} \mathord{\left/ {\vphantom {{(t)} {\partial {{p}_{i}}}}} \right. \kern-0em} {\partial {{p}_{i}}}}$ индекс l указывает векторную компоненту, индекс i – номер параметра, по которому выполняется дифференцирование. Приведем формулы для расчета псевдопроизводных в системе нормальных уравнений и расчета некоторых настоящих производных при этом (ср. выписанные выражения для ${{\partial {{q}_{i}}} \mathord{\left/ {\vphantom {{\partial {{q}_{i}}} {\partial {{p}_{j}}}}} \right. \kern-0em} {\partial {{p}_{j}}}}$ c кинематическими уравнениями в (4) и формулами (8))

Значения псевдопроизводных ${{\partial {{\varphi }_{l}}} \mathord{\left/ {\vphantom {{\partial {{\varphi }_{l}}} {\partial {{p}_{j}}}}} \right. \kern-0em} {\partial {{p}_{j}}}}$ $(j = 1,2, \ldots 6)$ определяются в процессе интегрирования уравнений

совместно с уравнениями (4), (5). Ненулевые начальные условия для ${{\partial {{\varphi }_{l}}} \mathord{\left/ {\vphantom {{\partial {{\varphi }_{l}}} {\partial {{p}_{j}}}}} \right. \kern-0em} {\partial {{p}_{j}}}}$ в точке ${{t}_{a}}$ равны 1 при ${{p}_{j}} = {{\theta }_{l}}$ и равны 0 в остальных случаях.

Прибавление найденных поправок $\Delta {{q}_{j}}({{t}_{a}})$ к имеющимся значениям ${{q}_{j}}({{t}_{a}})$ нарушает условие (7), поэтому новый кватернион ориентации нормируется. Внесенные нормировкой изменения уточненных компонент кватерниона являются величинами второго порядка относительно $\Delta {{q}_{j}}({{t}_{a}}).$

Точность аппроксимации измерений и оценки ${{z}_{ * }}$ будем характеризовать, следуя методу наименьших квадратов, соответствующими стандартными отклонениями. Последние вычисляются в предположении, что ошибки в измерениях некоррелированы и имеют одинаковые дисперсии, средние значения ошибок в измерениях одной и той же векторной компоненты напряженности МПЗ равны. Такой подход выбран из соображений удобства и вида функционала (6). При сделанных допущениях $z{\text{*}}$ – случайный вектор, который имеет приблизительно нормальное распределение со средним значением, равным истинному значению z. Вследствие условия (7) это распределение – несобственное, т.е. имеет вырожденную ковариационную матрицу. Чтобы избежать вырождения и сделать характеризацию ошибок более наглядной, ошибки $\Delta {{q}_{j}}({{t}_{a}})$ в задании компонент ${{q}_{j}}({{t}_{a}})$ вектора $z{\text{*}}$ представим в виде (8), где теперь ${{\theta }_{i}}$ образуют случайный вектор бесконечно малого поворота. Величины ${{\theta }_{i}}$ имеют нулевые математические ожидания и вместе с ошибками остальных компонент $z{\text{*}}$ описываются ковариационной матрицей

Здесь ${{\sigma }^{2}}$ – оценка дисперсии ошибок в измерениях (1) магнитометра 1, G – матрица $\left\| {{{G}_{{ij}}}} \right\|,$ вычисленная в точке $z{\text{*}}$. Точность аппроксимации измерений будем характеризовать полученным значением σ, точность оценки $z{\text{*}}$ – стандартными отклонениями $\sqrt {{{K}_{{11}}}} ,\sqrt {{{K}_{{22}}}} , \ldots \sqrt {{{K}_{{99}}}} .$ Стандартные отклонения величин ${{\theta }_{i}},$ ${{\chi }_{i}},$ ${{\gamma }_{1}},$ ${{\alpha }_{1}}$ и ${{\beta }_{1}}$ обозначим соответственно ${{\sigma }_{{\theta i}}},$ ${{\sigma }_{{\chi i}}},$ ${{\sigma }_{{\gamma 1}}},$ ${{\sigma }_{{\alpha 1}}}$ и ${{\sigma }_{{\beta 1}}}.$

Описанный подход практически без изменения применим в случае, когда вместо данных магнитометра 1 обрабатываются взвешенные данные (3) или данные магнитометра 2, пересчитанные в систему ${{y}_{1}}{{y}_{2}}{{y}_{3}}$ по формулам (2). Реализован также вариант обработки данных магнитометра 2 без пересчета их в систему ${{y}_{1}}{{y}_{2}}{{y}_{3}}.$ В этом варианте использовалась матрица перехода от системы ${{x}_{1}}{{x}_{2}}{{x}_{3}}$ к системе ${{z}_{1}}{{z}_{2}}{{z}_{3}}$ и уточнялись параметризирующие ее углы. Однако вариант обработки взвешенных данных дает наименьшее значение σ, и ниже приводятся результаты только этого варианта.

Вернемся к описанию методики. Чтобы обеспечить надежную сходимость описанного процесса, надо предусмотреть возможность его регуляризации. В данном случае регуляризация сводилась к предварительному использованию метода Левенберга–Марквардта [8] перед переходом к методу Гаусса–Ньютона. Почти всегда метод Левенберга–Марквардта плавно трансформировался в метод Гаусса–Ньютона, но иногда по окончании его работы метод Гаусса–Ньютона расходился. В таком случае в качестве $z{\text{*}}$ принимался результат, полученный методом Левенберга–Марквардта; он вполне обеспечивал требуемую точность реконструкции.

Пример реконструкции вращательного движения спутника на интервале времени 2 приведены на рис. 4–7. На рис. 4а изображены графики углов, характеризующих ориентацию спутника. Здесь λ и $\vartheta $ – углы отклонения оси ${{x}_{2}}$ от орта ${\mathbf{s}}$ направления “Земля–Солнце”, ${\mathbf{s}} = ({{S}_{1}},{{S}_{2}},{{S}_{3}})$ в системе ${{Y}_{1}}{{Y}_{2}}{{Y}_{3}}.$ Углы вычисляются по формулам

Плоскость отсчета угла $\vartheta $ содержит ось ${{Y}_{3}}$ и прямую “Земля–Солнце”. Направление отсчета этого угла – на север. Плоскость отсчета угла λ параллельна плоскости ${{Y}_{1}}{{Y}_{2}},$ направление отсчета угла λ – на восток. Угол φ – это угол между осью ${{x}_{1}}$ приборной системы и плоскостью орбиты; при $0 < \varphi < \pi $ эта ось лежит в полупространстве, в которое направлен кинетический момент орбитального движения спутника. Функция $\zeta (t)$ (нижний график) характеризует тень Земли: в тени $\zeta < 0,$ на свету $\zeta > 0.$ Эта функция рассчитывается по формуле

где ${\mathbf{r}}(t)$ – радиус-вектор центра масс спутника, ${{R}_{{\text{E}}}} = 6378.14$ км – радиус Земли, принятой шаром. Чтобы пояснить, на каких отрезках времени спутник находился в тени Земли, рядом с графиком функции $\zeta (t)$ проведена прямая $\zeta = 0.$ Как видно из графиков углов λ и $\vartheta $, спутник на большей части каждого обработанного интервала находился в режиме солнечной ориентации, совершая значительные развороты на коротких отрезках времени.

Рис. 4б содержат графики функций (5). Рис. 5а и 5б характеризуют построенную аппроксимацию магнитных измерений. На рис. 5а в каждой системе координат изображены две линии, которые практически сливаются. Одна из них – график функции ${{\hat {h}}_{i}}(t)$ (см. (6)), а другая – ломаная, звенья которой соединяют соседние по времени точки $({{t}_{n}},\tilde {g}_{i}^{{(n)}} - {{\Delta }_{i}}),$ $n = 1,2, \ldots N.$ Здесь $\tilde {g}_{i}^{{(n)}}$ – включенные в обработку измерения (3) при значениях ρ, ук-азанных в п. 2. На рис. 5б изображены аналогичные ломаные, соединяющие точки $({{t}_{n}},\Delta h_{l}^{{(n)}}),$ $\Delta h_{i}^{{(n)}} = \tilde {g}_{i}^{{(n)}} - {{\hat {h}}_{i}}({{t}_{n}}) - {{\Delta }_{i}}.$ Рис. 6 и 7 в увеличенном масштабе воспроизводят фрагменты графиков на рис. 4. Рис. 6 иллюстрирует развороты спутника, рис. 7 – режим солнечной ориентации. Судя по графикам углов λ и $\vartheta $, ориентация на Солнце оси ${{x}_{2}}$ спутника поддерживалась с требуемой точностью. Графики, иллюстрирующие реконструкции движения спутника на интервалах 1 и 3 приведены в [10].

Некоторые числовые характеристики результатов реконструкции на интервалах 1–3 в табл. 6. В табл. 7, 8 приведены оценки смещений в измерениях угловой скорости и оценки углов, задающих матрицу B (в скобках рядом с оценками указаны их стандартные отклонения).

Совместный анализ стандартных отклонений и собственных векторов матрицы G (см. п. 3), отвечающих ее нескольким минимальным собственным числам, позволяет понять характер неопределенности, которая возникает при реконструкции движения. Анализ основан на формуле

где ${{\mu }_{k}}$ и ${{u}_{k}}$ – собственные числа и ортонормированные собственные векторы матрицы G. Из этой формулы следует, что наибольшие стандартные отклонения имеют те определяемые параметры, которым отвечают наибольшие компоненты векторов ${{{{u}_{k}}} \mathord{\left/ {\vphantom {{{{u}_{k}}} {\sqrt {{{\mu }_{k}}} }}} \right. \kern-0em} {\sqrt {{{\mu }_{k}}} }}.$ Пусть собственные числа ${{\mu }_{k}}$ упорядочены по возрастанию: ${{\mu }_{1}} \leqslant {{\mu }_{2}} \leqslant \ldots \leqslant {{\mu }_{9}}$ и ${{\mu }_{1}} \ll {{\mu }_{2}}.$ В таком случае достаточно проанализировать только компоненты вектора ${{u}_{1}}.$ Для обработанных интервалов собственные числа ${{\mu }_{1}},$ ${{\mu }_{2}}$ и компоненты вектора ${{u}_{1}}$ приведены в табл. 9 (в расчетах использовались единицы измерения времени ${\text{1000}}$ с, напряженности магнитного поля 50 000γ; числа ${{\mu }_{7}},{{\mu }_{8}},{{\mu }_{9}}$ имеют порядок 5). Как видно из этой таблицы, ${{\mu }_{1}} \ll {{\mu }_{2}},$ компоненты ${{\chi }_{1}},{{\chi }_{2}},{{\chi }_{3}}$ весьма малы, а компоненты ${{\theta }_{1}},{{\theta }_{2}},{{\theta }_{3}}$ примерно совпадают с компонентами ${{\gamma }_{1}},{{\alpha }_{1}},{{\beta }_{1}}$. Эти соотношения поясняют малость стандартных отклонений ${{\sigma }_{{\chi {\kern 1pt} i}}}$ в табл. 7 и примерное совпадение стандартных отклонений ${{\sigma }_{{\theta 1}}}$ и ${{\sigma }_{{\gamma 1}}},$ ${{\sigma }_{{\theta 2}}}$ и ${{\sigma }_{{\alpha 1}}},$ ${{\sigma }_{{\theta 3}}}$ и ${{\sigma }_{{\beta 1}}}$ в табл. 6, 8. Такое совпадение говорит о связи ошибок в определении начальных условий ${{q}_{j}}({{t}_{a}})$ и углов ${{\gamma }_{1}},{{\alpha }_{1}},{{\beta }_{1}}.$

Обработка сокращенных интервалов измерений показала, что для надежного определения углов ${{\gamma }_{1}},$ ${{\alpha }_{1}},$ ${{\beta }_{1}}$ необходимо рассматривать интервалы, содержащие некоторое число интенсивных разворотов спутника.

Данная работа выполнена при поддержке РФФИ (проект 17-01-00143).