Содержание

1.     Введение………………………………………………………………………...4

2.     Шкалы времени…………………………………………………………………5

3.     Системы координат…………………………………………………………….6

4.     Модель движения Фобоса……………………………………………………...8

5.     Уравнения движения Фобоса………………………………………………….9

6.     Расчёт возмущающих ускорений…………………………………………….13

7.     Уравнения в вариациях……………………………………………………….18

8.     Значения астрономических постоянных…………………………………….22

9.     Техника определения параметров……………………………………………23

10.   Литература…………………………………………………………………….25

1.     Введение

На 2009 год в России запланирован полёт КА Фобос-Грунт с целью доставки на Землю образцов грунта спутника Марса Фобоса. Вблизи Фобоса КА предстоит выполнить сложные динамические операции [1-5]. При этом КА будет находиться в сфере действия гравитационных полей Солнца, Марса и Фобоса. Определение траектории КА проводится по совокупности траекторных измерений, выполненных с наземных станций слежения и измерений автономных систем на борту КА относительно Фобоса. В качестве наземных измерений используются радиотехнические измерения допплеровского смещения частоты  и времени распространения сигнала. С борта КА возможны оптические наблюдения Фобоса телевизионной системой и измерения расстояния до поверхности Фобоса  лазерным локатором. Поскольку орбита Фобоса известна недостаточно точно для решения задачи посадки, её определение проводится совместно с определением орбиты КА. Таким образом, для определения параметров движения КА относительно Фобоса необходимы алгоритмы для совместного уточнения начальных условий КА и Фобоса в рамках единой динамической модели движения. После выведения КА на орбиту спутника Марса, на орбите наблюдения Фобоса и на орбите вблизи Фобоса (квазисинхронная орбита) [2] процесс наведения аппаратуры КА на Фобос  и последующий процесс уточнения орбит требуют априорного знания о параметрах  орбиты Фобоса с достаточно высокой точностью. Такое знание должно быть получено заранее с использованием всех современных и исторических наблюдений Фобоса.

При построении теории использованы оптические измерения на интервале времени  1877-1989 годы , телевизионные измерения Фобоса с борта КA Mariner-9 (1971-1972 годы), Viking  (1976-1978 годы), Фобос-2 (1989 год), Mars Express (2005 год), измерения лазерного дальномера (MOLA) с борта Mars Global Surveyor (MGS) (1998 год), измерения углового расстояния между центром Фобоса и центром солнечного диска при наблюдениях с находящихся на поверхности Марса американских роверов Spirit и Opportunity (2004 год). Траекторные измерения КА, находящегося вблизи Фобоса, содержат важную информацию о его гравитационной постоянной. Такие измерения были использованы для КА Фобос-2 и MGS, поскольку они сближались с Фобосом и под воздействием его гравитации меняли траекторию. Виды измерений показаны на рис. 1.

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

2.     Шкалы времени

При построении теории используются две шкалы времени: TDB и UTC. Эфемериды Фобоса привязаны к шкале равномерного времени TDB, которое используется также в качестве аргумента дифференциальных уравнений для описания движения всех тел Солнечной систем. Привязка измерений из разных источников к моментам времени должна быть осуществлена к этой же шкале. Это также актуально для древних измерений, которые, как правило, были привязаны к времени, определяемом вращением Земли UT1. Разница между равномерным (эфемеридным) временем TDB и неравномерным временем UT1 определяется  эфемеридной поправкой . Для оптических измерений позапрошлого и прошлого века она приведена в работе Morley вместе со списком  источников измерений. Измерения последних десятилетий привязаны к гражданскому времени UTC, которое является приближением шкалы  времени UT1. Связь между координированным временем UTC и динамическим TDB задаётся линейным соотношением    в секундах,  где  (в юлианских датах) является -м моментом  изменения коэффициентов  или . До 1972 года применялось лусочно-линенйное представление, теперь все , а  изменяются на 1 сек, когда разность UT1-UTC достигает 0.9 . Это делается либо в начале, либо в середине года.

Юлианский счет суток используются для подсчета числа суток, протекших между заданными датами и других расчетов, связанных с привязкой событий ко времени. Юлианские даты считаются непрерывно через годы, столетия и тысячелетия, начиная с 12 часов гринвичского времени 1 января 4713 г. до нашей эры. Этот момент времени считается нулевым юлианским днем. Существуют алгоритмы для вычисления номера юлианского дня по календарной дате и обратно.

3.     Системы координат

В работе для описания движения небесных тел и Фобоса используются системы координат (СК) двух видов: инерциальная СК и СК, связанные с фигурой небесного тела. Для описания движения Фобоса применяется ареоцентрическая геоэкваториальная СК, связанная с экватором и равноденствием Земли эпохи J2000. Далее для обозначения этой системы координат используется аббревиатура IMCS (Inertial Mars centric coordinate system) или просто J2000. Необходимо переходить из СК, жестко связанных с фигурами Земли  (гринвическая система координат GCS) , Марса (Magr), и Фобоса (Phgr) в инерциальную СК и обратно.

Для определения СК, связанных с фигурами Марса и Фобоса используются угловые координаты северного полюса и положение первого меридиана. Угловые координаты северного полюса – прямое восхождение  и склонение  – отсчитываются от осей фундаментальной системы координат J2000, Положение первого меридиана определяется углом , который отсчитывается вдоль экватора планеты в восточном направлении от точки пересечения плоскости экватора Земли СК J2000 с плоскостью экватора Марса или Фобоса до точки пересечения первого меридиана с экватором.

Переход от СК J2000 к GCS

Если положение объекта в СК J2000 определено вектором , то его положение в GCS выражается формулой , где матрица  перехода из СК J2000 в GCS  вычисляется как произведение четырех матриц ,

 – матрица прецессии, определяющая переход от среднего экватора и равноденствия фундаментальной эпохи J2000 к среднему экватору и равноденствию текущей эпохи ;

 – матрица нутации в эпоху ;

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

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

Обратное преобразование выполняется по формуле, где

 – матрица, транспонированная по отношению к . Параметры для вычисления матриц определены решениями Международного астрономического союза (МАС) [6].

Переход от J2000 к СК, связанной с небесным телом (Марсом и Фобосом)

Пусть  обозначает координаты объекта в геоэкваториальной геоцентрической СК J2000, тогда его положение  (CBCS – Celestial Body Coordinate System) в СК небесного тела (НТ), связанной с экватором и нулевым меридианом НТ определяются формулой , где  – положение центра масс НТ относительно центра масс Земли в СК J2000,  – матрица перехода от СК J2000 к СК Марса или Фобоса.

Если в момент времени  известны координаты северного полюса НТ  и угол , то матрицу перехода  можно представить как произведение двух матриц , где

Формулы для вычисления , ,  предписаны решениями МАС  [6].

4.     Модель движения Фобоса

Модель движения Фобоса учитывает возмущающие ускорения, вызванные нецентральностью поля тяготения Марса, гравитационное влияние Солнца, Деймоса, планет и Луны, вековое ускорение Фобоса, возникающее из-за воздействия приливных сил, а также влияние гравитационного поля Фобоса на его поступательное движение вокруг Марса. Вековое ускорение в движении Фобоса по орбите обнаружил в 1945 году американский астроном Б. Шарплесс. Спутник настолько близок к Марсу, что его перемещение вызывает приливы в марсианской коре. Приливная волна из-за трения и сил напряжения в коре планеты отстает от движения Фобоса примерно на четверть круга. Притяжение отстающих приливных горбов тормозит движение спутника, в результате чего Фобос движется по очень пологой спирали, постепенно приближаясь к поверхности Марса. Время баллистического существования Фобоса составляет приблизительно 15 миллионов лет [7].

Расчет кинематических параметров Фобоса осуществляется численным интегрированием уравнений его движения в элементах, чем достигается увеличение скорости интегрирования по сравнению с интегрированием в прямоугольных координатах. Уравнения учитывают все упомянутые выше возмущения, за исключением воздействия приливных сил. Традиционно [8,9] учет векового ускорения, вызванного марсианскими приливами, осуществляется путем изменения долготы Фобоса. В настоящей работе вековое смещение Фобоса вдоль своей траектории производится изменением момента времени  прохождения ближайшего к текущему моменту времени  узла орбиты, на величину . Такой подход полностью эквивалентен традиционному изменению долготы, но больше соответствует использующейся форме записи уравнений движения и набору орбитальных элементов. Коэффициент  связан с вековым ускорением в долготе формулой .

5.     Уравнения движения Фобоса

В текущий момент времени  параметры движения Фобоса определяется либо его вектором состояния , либо кеплеровскими элементами орбиты , где

 – координаты центра масс Фобоса в IMCS,

 – компоненты вектора скорости в IMCS,

 _– большая полуось,

,  компоненты вектора Лапласа,

– наклонение,

 – долгота восходящего узла,

 – время прохождения восходящего узла,

угловые параметры  отсчитываются от осей IMCS.

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

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

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

Частные производные пертурбационной функции  по оскулирующим элементам орбиты  вычисляются по формулам

,

где

 – вектор (столбец) возмущающего ускорения движения Фобоса,

 – вектор, транспонированный по отношению к

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

1.     По элементам орбиты  вычисляется вектор состояния  для этого

1.1.       Определяются

среднее движение ,

средняя аномалия ,

эксцентриситет ,

параметр орбиты ,

где

 – текущий момент времени

 – гравитационная постоянная Марса,

 – гравитационная постоянная Фобоса.

1.2.       Определяется эксцентрическая аномалия  как решение уравнение Кеплера

Решение проводится итерационно методом Ньютона

 (Начальное приближение)

1.3.       Вычисляются синус и косинус истинной аномалии

,  и

модуль радиуса-вектора

1.4.       Вычисляются орбитальные координаты

, ,

и скорости

, ,

1.5.       Формируется матрица  перехода от орбитальных координат  к координатам IMCS

1.6.       Вычисляется вектор состояния Фобоса

,

2.      Элементы матрицы  вычисляются по формулам

где

 – единичные векторы орбитальной системы координат

 – аргумент широты,

 – истинная аномалия.

6.     Расчет возмущающих ускорений

Возмущающее ускорение Фобоса  в системе IMCS складывается из ускорения, вызванного нецентральностью поля тяготения Марса , возмущений от Солнца и планет солнечной системы , Деймоса , ускорения, вызванного влиянием поля Фобоса на его поступательное движение:

.

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

Возмущающая функция, обусловленная нецентральностью гравитационного поля, в системе MCS имеет вид:

.где

– радиус-вектор точки с прямоугольными координатами ,

 – долгота и широта этой точки,

 – гравитационная постоянная Марса

 – экваториальный радиус Марса,

 – коэффициенты при зональных гармониках,

 – коэффициенты при тессеральных и секториальных гармониках,

 – полиномы Лежандра,

 – присоединённые функции Лежандра.

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

,

где  – m-ая производная полинома Лежандра ,

.

Присоединенные функции Лежандра представляются формулой

      или     .

При этом  ,  ,  .

Ускорение, обусловленное нецентральностью гравитационного поля, равно градиенту возмущающей функции :

 (в прямоугольной системе координат).

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

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

, где  – символ Кронекера.

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

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

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

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

Здесь

 гравитационная постоянная тела

 координаты Фобоса в эпохе J2000

 координаты возмущающего тела  в эпохе J2000.

Для вычисления координат Солнца и планет   используются таблицы эфемерид DE405 .

Возмущения от второго спутника Марса Деймоса вычисляются по той же формуле, но вследствие незначительности этих возмущений, для получения координат Деймоса на интервале с 1870 по 2030 (для насчитывания эфемериды) годы  используется  его кеплерова орбита, которая периодически обновляется с заданным шагом. Для этого используются данные теории Chapront-Touzé [9]. Ниже даётся таблица гравитационных постоянных тел, возмущения от которых учитываются в модели. Эти величины соответствуют значениям, принятым в теории DE405, за исключением Марса, значение которого согласовано с коэффициентами гравитационного поля. Размерность .

Вычисление   происходит следующим образом. В фобосоцентрической инерциальной системе J2000 определяются координаты Марса . Далее они переводятся в фобосографическую систему координат . Описание матрицы перехода  приводится в разделе 3. По вышеописанной методике вычисляется ускорение , вызванное воздействием гравитационного поля Фобоса на Марс. Используются коэффициенты, приведённые в разделе 8. Ускорение в точке, совпадающей с центром Фобоса, вычисляется по формуле

где  гравитационная постоянная Фобоса.

7. Уравнения в вариациях

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

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

,

,

,

,

.

Здесь , .

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

При вычислении вклада от Марса введём единичный вектор . Вычисляется множитель  и далее матрица

Выполняется умножение  на .

,

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

,

Эта часть возмущений от Солнца и планет проводилась в эпохе J2000.

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

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

, .

Далее необходимо получить

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

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

8.     Значения астрономических постоянных

В настоящей работе использовано поле Марса 8x8 [10], и нижеследующая таблица отражает величины принятых констант. Экваториальный радиус Марса равен 3397.0 км

Таблица ненормализованных коэффициентов гравитационного поля Марса

Таблица гравитационных постоянных тел в модели движения Фобоса

Гравитационное поле Фобоса представлено только сжатием, величина которого  взята из работы [11] пересчитана в ненормализованный вид:

С20= - 0.104, экваториальный радиус, используемый в расчётах 12 км.

8.     Техника определения параметров

При построении теории наряду с элементами Фобоса уточнялись параметры движения тех КА, которые близко подходили к Фобосу, - Фобоса-2 и MGS. В число уточняемых параметров   включены: элементы орбиты Фобоса , гравитационная постоянная Фобоса , вековое ускорение Фобоса , параметры движения КА Фобос-2  и параметры манёвров, параметры движения КА MGS .

Определение параметров теории  методом максимального правдоподобия [12] требует минимизации функционала где , и  является  м вектором-измерением с ковариационной матрицей априорных ошибок  и ,  является вычисляемым значением вектора-измерения. Для минимизации  применяется итерационный метод. На каждом - ом шаге решается линейная система

 и параметры  уточняются.

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

1.     Подсистема Фобос, ;

2.     Подсистема Фобос-2 (с уточнением манёвров), Фобос, ;

3.     Подсистема MGS, Фобос

Каждая -я подсистема представляет собой задачу обработки определённого набора измерений на каждой итерации. Она начинается подготовкой необходимых данных – измерений, начальных условий и др. и заканчивается формированием матрицы  и правых частей  нормальных уравнений для этой подсистемы , которые войдут в общую систему в качестве составных частей. Первая подсистема обеспечивает обработку всех наземных оптических измерений Фобоса и прохождения Фобоса по диску Солнца. Она производит расчёт параметров движения Фобоса на интервале более 100 лет и обновляет начальные условия для двух других подсистем перед началом каждой следующей итерации. Вторая подсистема проводит обработку траекторных измерений радиальной скорости КА Фобос-2 и бортовых телевизионных измерений Фобоса. Область действия подсистемы распространяется на интервал от начала переходной орбиты, на которую КА вышел после проведения манёвра торможения и перехода на орбиту спутника Марса, и закачивается квазисинхронной орбитой (конец активного существования аппарата). Описание движения КА на всех участках сделано в рамках одной траектории, на которой есть участки пассивного движения и маневрирование. Эту траекторию описывают элементы орбиты КА в начальный момент времени на переходной орбите и параметры манёвров, а также значение гравитационной постоянной Фобоса и его вековое ускорение. Каждый манёвр описывается шестью параметрами, обеспечивающими расчёт приращений вектора состояния  при завершении манёвра. В обработку включены априорные значения параметров манёвров с учётом ошибок их исполнения. Эта априорная информация представлена в форме ковариационной матрицы. Третья подсистема предназначена для обработки измерений MGS. Протяжённость интервала измерений взятого в обработку составляет несколько суток, модель движения описывается в рамках пассивного движения.

Схематически  формирование общей системы представлено на  рис. 2. Каждый элемент общей системы получается суммированием соответствующих элементов подсистем, что показано одинаковым окрасом. Кроме того, на рис. 2 показано место будущей подсистемы КА3 Фобос-Грунт и параметры для уточнения из этой подсистемы.

Литература

1.     Аким Э.Л., Заславский Г.С.,Морской И.M. и др. Баллистика, навигация и управление полётом космического аппарата в проекте «Фобос-Грунт» // Изв. РАН ТиСУ. 2002. 5. с. 153-161

2.     Тучин А.Г. Квазисинхронные орбиты и их использование для сближения космического аппарата с Фобосом // Космические исследования 2007. 45. 2. с. 144- 149

3.     Akim E.L., Botkin A.V., Stepaniants V.A. et al. Orbit selection, Navigation and Maneuvers before the Landing on the Phobos Surface for Phobos Sample Returne Project // The Preceedings of the 17^th International Symposium on Soace Flight Dynamics, 16-20 June 2003. Moscow. Russia. V.1.

4.     Akim E.L. Popove G.A. Tuchin A.G. Mechanics and Motion Control of a Space Vehicle in the Project of Relict Substance Delivery on Earth (The Project “Phobos-Grunt”) // 16^th IFAC Symposium on automatic control in Aerospace. Preprints, Saint Peterburg. Russia. 2004. V.1

5.     AkimE.L., Stepaniants V.A., Shishov V.A. et al. // Ballistics,navigation and motion control of the SC on stages of the Phobos surface approaching and landing // The proceedings of the 18^th International Symposium on Space Flight Dynamics, 11-15 October 2004. Munich, Germany

6.     Seidelmann P.K., Abalakin V. K., Bursa M. et al. Report of the IAU/IAG working group on cartographic coordinates and rotational elements of the planets and satellites:2000 // Celestial Mechanics and Dynamical Astronomy. 2002. 82. P. 83-110

7.     Veverka J., Burns J. A. The moons of Mars // Annual review of earth and planetary sciences 1980. 8. P. 527-558

8.     Chapront-Touze M. ESAPHO: a semi-analytical theory for the orbital motion of Phobos // Astronomy and Astrophysics. 1988. 200. P. 255-268

9.     Chapront-Touze M. Orbits of the Martian satellites from ESAPHO and ESADE theories // Astronomy and Astrophysics. 1990. 240. P. 159-172

10.     F. G. Lemoine, D. E. Smith, D.D. Rowlands, M.T. Zuber, G. A. Neumann, and D. S. Chinn, An improved solution of the gravity field of Mars (GMM-2B) from Mars Global Surveyor, J. Geophys. Res., 106(E10), 23359-23376, October 25, 2001.

11.     Martinec Z., Pec K. The Phobos gravitational field modeled on the basis of its topography // Earth, Moon, and Planets. 1989. 45. P. 219-235

12.      Аким Э.Л., Энеев Т.M. Определение параметров движения космического аппарата по данным траекторных измерений // Космические исследования. 1963. 1. Вып. 1. с. 5-50

Рис.1 Типы измерений, использованных для построения теории.

Рис. 2 Схема формирования общей системы нормальных уравнений  для получения параметров  на очередной итерации.