Определение параметров вращательного движения КА «Монитор-Э» по телеметрическим данным о токе солнечных батарей

( Reconstruction of the spacecraft Monitor-E attitude motion by measurements of the current from solar arrays
Preprint, Inst. Appl. Math., the Russian Academy of Science)

1. Математическая модель вращательного движения КА. КА считаем твёрдым телом, геоцентрическое движение центра масс которого – кеплерово эллиптическое. Элементы этого движения находятся по данным радиоконтроля орбиты. Для записи уравнений движения КА относительно центра масс введём две правые декартовы системы координат – орбитальную  и образованную главными центральными осями инерции КА . Точка  – центр масс КА, оси  и  направлены соответственно вдоль геоцентрического радиус-вектора точки  и по трансверсали к орбите в этой точке.

         Положение системы  относительно системы  зададим углами ,  и , которые введём следующим образом. Система  может быть переведена в систему  тремя последовательными поворотами: 1) на угол  вокруг оси , 2) на угол  вокруг новой оси , 3) на угол  вокруг новой оси , совпадающей с осью . Матрицу перехода от системы  к системе  обозначим , где – косинус угла между осями  и . Элементы этой матрицы выражаются через введённые углы с помощью формул

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

Здесь  и  – соответственно компоненты векторов абсолютной угловой скорости КА и скорости его центра масс относительно поверхности Земли в системе ;  – среднее орбитальное движение КА; – моменты инерции КА относительно осей ,  – плотность атмосферы в точке ;  – параметры аэродинамического момента. Плотность атмосферы рассчитывается согласно модели [1].

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

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

Ток, вырабатываемый СБ, примерно пропорционален косинусу угла падения солнечных лучей на их светочувствительную поверхность. Последняя расположена в плоскости, неподвижной относительно системы . Нормаль  к светочувствительной стороне СБ задается двумя углами: углом  между вектором  и его проекцией на плоскость   и углом  между осью  и проекцией  на плоскость  . Направление отсчёта  согласовано с направлением оси . Компоненты  в системе  вычисляются по формулам: , , . Орт  направления «Земля-Солнце» задается в орбитальной системе координат компонентами , которые рассчитываются по приближенным формулам [2] и элементам кеплеровой орбиты КА. Упомянутый косинус и снимаемый с СБ ток задаются формулами

 Здесь  – ток, вырабатываемый СБ на орбите Земли при перпендикулярном падении солнечных лучей на их плоскость, А.

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

         Телеметрическая информация о токе, снимаемом с батарей, представляет собой последовательность чисел

Здесь  – приближённое значение тока в момент времени , . Разности , как правило, не превышают нескольких секунд. В обработку включаются отрезки данных, длина которых  составляет от 20 до 40 минут.

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

Аппроксимацией фактического движения КА на этом отрезке будем считать решение, доставляющее такому функционалу минимум. Минимизация  проводится по начальным значениям параметров движения КА в точке : , , ,  и параметрам , , , , . Для простоты письма объединим все уточняемые величины в один вектор  и будем рассматривать функционал (4) как функцию . Тогда  – искомая оценка вектора .

         Минимизация функционала (4) выполнялась в два этапа: сначала методом случайного поиска находилось грубое приближение , которое затем уточнялось методом Левенберга-Марквардта.

3. Реализация метода Левенберга-Марквардта. Применение этого метода в задачах определения вращательного движения спутников по данным измерений бортовых датчиков описано в [3, 4]. Метод является одним из вариантов метода Гаусса-Ньютона. На каждой итерации этого метода поправка , уточняющая имеющееся приближённое значение , определяется системой

Здесь  – матрица нормальных уравнений, – положительный параметр,   – единичная матрица порядка 13,  – символ Леви-Чивиты,  – псевдопроизводные [3], служащие для представления истинных производных  при . Значения   определяются в процессе интегрирования уравнений

совместно с уравнениями (1) и уравнениями в вариациях относительно . Ненулевые начальные условия  и  имеют вид

,    ,    ,    .

Ненулевые производные  в уравнениях (5) определяются следующим формулами:

,       ,

,     ,     .

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

,

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

         Чтобы в результате минимизации функционала (4) значения параметров , ,  лежали в разумных с физической точки зрения пределах, в этот функционал вводились дополнительные слагаемые

.

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

4. Поиск начального приближения. Чтобы алгоритм минимизации Левенберга-Марквардта был надёжен, необходимо иметь достаточно точное начальное приближение точки минимума. Для этого использовался метод случайного поиска с обучением [5].

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

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

Определить значение  для КА, повёрнутого на произвольные углы  (вокруг орта ) и затем  (вокруг орта ), можно по формулам

, , ,

Здесь матрица  определяется положением Солнца в орбитальной системе координат в начальный момент времени, матрица  определяется положением  в системе координат  (раздел 2),  – матрица последовательных поворотов КА: вокруг орта  на угол  и затем вокруг орта  на угол . Элементы матрицы  зависят также от начального угла  между ортами  и .

Соотношение (6) позволяет выразить углы ,  и  в функции углов  и . Эти углы ( и ) использовались в процедуре случайного поиска для задания начальной ориентации КА. Кроме того, для удобства выбора начального значения и последующего контроля величины угловой скорости КА, компоненты  () вектора  были заменены величинами  , , .

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

Алгоритм случайного поиска с обучением [5] состоит в построении последовательности векторов  и  по следующему правилу. На первом шаге в качестве  берётся произвольно, и принимается . На каждом последующем шаге вычисляются векторы

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

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

После окончания случайного поиска по найденным компонентам вектора  определяется начальное приближение вектора  для вычислений методом Левенберга-Марквардта.

5. Результаты определения вращательного движения КА. Определение фактического движения КА относительно центра масс по данным (3) было выполнено на 12 интервалах времени. Основные характеристики этих интервалов приведены в табл. 1. Здесь для каждого интервала указаны: дата и декретное московское время первого измерения , длина интервала , число  включённых в обработку измерений. Полученные результаты представлены в табл. 2, 3 и на рис. 2¸15. В таблицах приведены результаты минимизации функционала (4) на интервалах из табл. 1. Здесь указаны значения некоторых компонент вектора , стандартные отклонения этих компонент и стандартное отклонение  ошибок данных (4). Параметры , ,  и  были определены со значительной погрешностью, поэтому в таблицах не приводятся. Погрешность определения начального углового положения КА на интервалах  характеризуется значениями , на интервалах 7¸12 – значениями , причем доминируют значения, близкие верхним границам. Стандартные отклонения параметров  лежала в пределах от 0.01 до 1.5. Большая погрешность определения угловых параметров , ,  вызвана, по-видимому, неблагоприятным расположением векторов  и , а параметров  – слабым влиянием аэродинамического момента на вращательное движение КА.

Анализ графиков найденных угловых скоростей КА показывает, что вращательное движение КА на интервалах 1¸6 близко к регулярной прецессии Эйлера с преимущественным вращательным движением КА вокруг оси . На интервалах 7¸12 характер движения КА более сложный, однако преимущественное вращение спутник также совершает вокруг оси .

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

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

Аналогичным образом углы  и  задают в системе  положение вектора  (рис. 2). Угол между векторами  и  обозначим .

Анализ результатов определения движения КА показал, что на интервале 1   и , на интервалах 2¸6   и . Эти величины существенно меньше угловой скорости движения КА вокруг , равной на интервалах  примерно °/с. На интервалах 7¸12   и , а угловая скорость КА равна °/с.

Для наглядного представления движения КА вокруг вектора , будем рисовать проекции годографа орта оси  – оси собственного вращения КА – на координатные плоскости системы , связанной с этим вектором [6]. Система  определяется следующим образом (рис. 3). Ось  направлена вдоль ; ось  лежит в плоскости , перпендикулярна  и образует тупой угол с осью ; ось  дополняет систему до правой. Помимо проекций годографа орта оси  движение КА будем характеризовать углом нутации , который эта ось образует с вектором .

Учитывая характер движения КА на рассматриваемых интервалах, это движение можно характеризовать усреднёнными величинами

,      .

При отсутствии внешних механических моментов угловые скорости КА-твёрдого тела относительно главных центральных осей инерции будут определяться по соотношениям , , , где ,  – постоянные величины,

,      .

Средние квадратические отклонения

,     

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

В табл. 4 приведены некоторые параметры вращательного движения КА для рассмотренных интервалов. Вращательное движение КА характеризуется здесь параметрами , ,  и , а также максимальным значением угла  на интервале – параметром . Движение векторов  и  происходило медленно и характеризуется усреднёнными значениями , , , , . Усреднение выполнялось по формулам, приведённым для параметров  и .

На рис. 4 – 15 представлена реконструкция вращательного движения КА на этих интервалах. Эти рисунки иллюстрируют решения уравнений вращательного движения КА, доставляющие минимум функционалу (4). Все рисунки скомпонованы одинаково и содержат проекции годографа орта  на плоскости системы , графики углов  и , графики угловых скоростей КА , ,  и графики расчетного тока СБ . Маркерами рядом с графиками  указаны данные измерений . Все графики построены на соответствующих отрезках ; началом отсчета времени служит точка . Угловые величины выражены в градусах, угловые скорости – в град./с, ток – в амперах.

         Исследуемые интервалы движения КА, можно разделить на интервалы быстрого движения КА (интервалы 1¸6, рис. 4 – 9) и медленного движения (интервалы 7¸12, рис. 10 – 15). На быстрых интервалах движение КА представляет собой регулярную прецессию Эйлера с преимущественным вращением вокруг оси . Кинетический момент КА на разных интервалах находился как вблизи орбитальной плоскости, так и далеко от последней. Скорость движения вектора кинетического момента была значительно меньше скорости движения КА. На более медленных интервалах движение КА носило более сложный характер, хотя проекции оси  на направление вектора  сохраняли свой знак. Как и следовало ожидать, в случае более медленного вращательного движения КА и более сложного изменения во времени тока, снимаемого с СБ, движение определялось лучше, чем в случае быстрого вращения КА вокруг оси . Эта ситуация является общей, хотя возможны исключения. Примером может служить случай, когда коллинеарны или достаточно близки по направлению вектор нормали к СБ, ось вращения КА и вектор, направленный на Солнце.

Данная работа выполнена в рамках проекта РФФИ 08-01-00467.

Литература.

Таблица 1. Интервалы определения вращательного движения КА.

Таблица 2. Результаты минимизации.

Таблица 3. Результаты минимизации.

Таблица 4. Параметры движения КА.