«ЧАЗОВ Вадим Викторович РАЗРАБОТКА И ПРИМЕНЕНИЕ АЛГОРИТМОВ ЧИСЛЕННО-АНАЛИТИЧЕСКОГО МЕТОДА ВЫЧИСЛЕНИЯ ПОЛОЖЕНИЙ ИСКУССТВЕННЫХ СПУТНИКОВ ЗЕМЛИ Диссертация на соискание учёной степени доктора физико-математических наук ...»

МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

имени М. В. ЛОМОНОСОВА

ГОСУДАРСТВЕННЫЙ АСТРОНОМИЧЕСКИЙ ИНСТИТУТ

имени П. К. ШТЕРНБЕРГА

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

УДК 521.13

ЧАЗОВ Вадим Викторович

РАЗРАБОТКА И ПРИМЕНЕНИЕ АЛГОРИТМОВ

ЧИСЛЕННО-АНАЛИТИЧЕСКОГО МЕТОДА

ВЫЧИСЛЕНИЯ ПОЛОЖЕНИЙ

ИСКУССТВЕННЫХ СПУТНИКОВ ЗЕМЛИ

Диссертация на соискание учёной степени доктора физико-математических наук

Специальность 01.03.01. Астрометрия и небесная механика Москва – 2012 Содержание Содержание Предисловие 1 Постановка задачи 1.1 Стандартные соглашения...................... 1.1.1 Системы отсчёта....................... 1.1.2 Пространство-время..................... 1.1.3 Окрестности Земли...................... 1.1.4 Солнечная система...................... 1.1.5 Прецессия и нутация..................... 1.1.6 Стандартное преобразование................ 1.1.7 Координаты обсерваторий.................. 1.2 Исходные данные........................... 1.2.1 Эволюция элементов орбиты................ 1.2.2 Уравнения движения..................... 1.2.3 Геопотенциал......................... 1.2.4 Численная модель...................... 1.2.5 Приливный потенциал.................... 1.3 “Стационарный” объект....................... 1.4 Промежуточный потенциал..................... 1.5 Канонические преобразования................... 1.6 “Осреднённые” уравнения...................... 1.7 Постановка задачи.......................... 2 Основные алгоритмы 2.1 Промежуточная орбита....................... 2.1.1 От вектора состояния к параметрам орбиты....... 2.1.2 Постоянные интегрирования 1, 2, 3......... Позиционные параметры a, e,, 2.1.3........... 2.1.4 Алгоритм операций с полиномами............. 2.1.5 Угловые переменные l0, g0, h0............... Содержание 2.1.6 От параметров орбиты к вектору состояния........ 2.1.7 Проверка алгоритмов.................... 2.1.8 Новые результаты...................... 2.2 Частные производные........................ 2.2.1 Переменные действия.................... 2.2.2 Частные производные по элементам L, G, H...... 2.2.3 Частные производные по параметрам a, e,...... 2.2.4 Производные от угловых переменных........... 2.2.5 Новые результаты...................... 2.3 Алгоритм преобразования...................... 2.3.1 Аномальный геопотенциал................. 2.3.2 Притяжение светил...................... 2.3.3 Эффект приливов...................... 2.3.4 Световое давление...................... 2.3.5 Положения светил...................... 2.3.6 Начальные функции..................... 2.3.7 Элементарное слагаемое................... 2.3.8 Сумма слагаемых....................... 2.3.9 Формулы рекурсии...................... 2.3.10 Новые результаты...................... 2.4 Алгоритм интегрирования...................... 2.4.1 Дифференциальные соотношения............. 2.4.2 Рекуррентный алгоритм................... 2.4.3 Особые слагаемые...................... 2.4.4 Новые результаты...................... 2.5 Осреднённые уравнения....................... 2.5.1 “Несингулярные” переменные................ 2.5.2 Частные производные.................... 2.5.3 Короткопериодические неравенства............ 2.5.4 Порядок действий...................... 2.5.5 Функции координат..................... Содержание 2.6 Дополнительные алгоритмы..................... 3.2 Вычислительные аспекты...................... 3.4 Сравнительные испытания..................... 3.5.4 Сравнение результатов, пункт 7090............. 3.5.5 Сравнение результатов, пункт 7825............. 3.5.6 Сравнение результатов, пункт 7840............. Простое дифференцирование....................... Полиномы и присоединённые функции Лежандра........... Лагеос и Лагеос-2. Состав наблюдений................. Лагеос и Лагеос-2. Поле скоростей.................... Список иллюстраций 1 Стелла, средние элементы орбиты................. 6 Лагеос и Лагеос-2, стандартное отклонение............ 7 Лагеос и IERS, разности координат полюса............ 8 Лагеос-2 и IERS, разности координат полюса........... 9 Данные Международной службы вращения Земли........ 10 Результаты обработки наблюдений Лагеос и Лагеос-2...... 11 Коэффициент эмпирического ускорения.............. 12 Эмпирический коэффициент отражения.............. 13 Лагеос, пункт 7090, “невязки”.................... 14 Лагеос, пункт 7090, “невязки”, линейная аппроксимация..... 15 Лагеос, пункт 7825, “невязки”.................... 16 Лагеос, пункт 7825, “невязки”, линейная аппроксимация..... 17 Лагеос, пункт 7840, “невязки”.................... 18 Лагеос, пункт 7840, “невязки”, линейная аппроксимация..... 19 Расчёт ситуаций сближения..................... 20 “Электро-Л”, прогноз ситуаций сближения на 2007 год...... 24 Молния 3-39, активный режим полёта............... 25 Молния 3-39, коррекция орбиты.................. 26 Молния 3-39, пассивный режим полёта.............. 27 Молния 3-39, фильтрация набюдений............... 29 Молния 3-39, данные “NORAD”................... 30 Глонасс 118, фильтрация наблюдений............... 31 Метеор-3М, фильтрация наблюдений............... 32 Блиц, фильтрация наблюдений................... 33 Ларец, фильтрация наблюдений.................. 34 Небесный промежуточный экватор................. Список таблиц 2 Числовые значения коэффициентов зональных гармоник.... 3 Коэффициенты многочленов от двух переменных........ 4 Функция r0 / как сумма элементарных слагаемых....... 5 Функция как сумма элементарных слагаемых......... 6 Начальные функции: число слагаемых............... 7 Функции координат: число слагаемых............... 8 Оценка погрешности прогноза на основе “средних” элементов.. 9 Производящая функция для полиномов Лежандра........ 10 Величины коэффициентов полиномов............... 14 Лагеос, сравнительные испытания, пункты 7080 и 7105..... 15 Лагеос, сравнительные испытания, пункт 7090.......... 16 Лагеос, сравнительные испытания, пункт 7941.......... 17 Лагеос, сравнительные испытания, несколько пунктов...... 18 Лагеос-2, сравнительные испытания................ 19 Эталон-1, сравнительные испытания................ 20 Обсерватория Яррагади, сравнительные испытания....... 21 Поле скоростей для пяти пунктов................. 22 Параметры орбиты объекта 95334................. 23 Глонасс 109, фильтрация наблюдений............... 24 Разность нуль-пунктов шкал времени............... 26 Функция 27 Функция 28 Функции координат: несколько слагаемых............. Предисловие В диссертации представлена новая численно-аналитическая теория движения искусственных спутников Земли (ИСЗ). Теория была применена для определения параметров движения космических объектов и получения достоверных оценок геодинамических параметров на основе наблюдений.

Численно-аналитический подход позволяет использовать преимущества как аналитических методов, так и метода численного интегрирования уравнений движения.

Классическая форма аналитического способа вычисления положений небесных объектов заключается в следующей процедуре: на вековые изменения параметров орбиты накладываются долгопериодические неравенства и короткопериодические возмущения. Преимуществом аналитического метода является возможность увеличения скорости расчёта положений объекта при условии существенных ограничений на точность прогноза.

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

Численно-аналитический подход объединяет оба метода: часть неравенств, имеющих “короткий” период, определяется аналитически, а долгопериодические, резонансные и вековые слагаемые возмущающей функции составляют эволюционный гамильтониан – основу численного интегрирования “осреднённых” уравнений движения. Шаг интегрирования “осреднённых” уравнений значительно превышает значение шага интегрирования, задаваемого в численных методах. Конкретное число, разделяющее слагаемые “короткого” и “долгого” периодов, является свободным параметром.

Запуск Первого искусственного спутника Земли вызвал интерес учёных к задаче определения орбит объектов, возмущаемых аномалиями геопотенциала, притяжением Луны и Солнца, потенциалом, обусловленным приливами упругой Земли, сопротивлением атмосферы и световым давлением.

Профессором Е.П.Аксёновым была построена аналитическая теория движения ИСЗ на основе решения обобщённой задачи двух неподвижных центров. Результаты представлены в монографии [5].

Было открыто и получило признание обобщение одного из методов теории возмущений – метода канонических преобразований [156].

Тогда же был успешно применён численно-аналитический метод расчёта спутниковых орбит [36]. Важные результаты в этом направлении получены профессором М.Л.Лидовым и его учениками, некоторые итоги подведены в обзорном докладе [96], там же намечены перспективы исследований.

Тем не менее, в наши дни все центры обработки высокоточных наблюдений ИСЗ проводят расчёты с помощью программ, в которых прогноз положений объектов выполняется алгоритмом численного интегрирования дифференциальных уравнений движения [109, 204].

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

В качестве исходных материалов предлагаемого исследования были использованы:

рекомендации Международного астрономического союза [199], стандартные соглашения Международной службы вращения Земли [189], лазерные наблюдения искусственных спутников Земли [190], позиционные наблюдения космических объектов [178], база данных об элементах орбит более 14000 объектов [206], база данных о точных положениях навигационных спутников [197].

Целью исследований является решение следующих задач:

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

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

дифференцирование и интегрирование слагаемых возмущающей функции;

составление и численное интегрирование “осреднённых” уравнений движения и учёт короткопериодических неравенств с использованием “несингулярных” элементов орбиты;

обработка высокоточных лазерных измерений топоцентрических дальностей до спутников Лагеос и Лагеос-2 на длительных интервалах времени и оценка значений геодинамических параметров;

разработка методики предсказания ситуаций сближения ИСЗ;

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

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

В тексте диссертации термин алгоритм объединяет несколько понятий:

это и связанные между собой формулы и соотношения, и последовательность действий, и реализация процедуры вычислений на компьютере. Совокупность алгоритмов, предназначенная для решения поставленной задачи, в тексте называется “программным обеспечением” или “пакетом программ”.

Для отладки программного обеспечения и оценки точностных характеристик алгоритмов была разработана методика сравнительных испытаний.

Методика основана на использовании высокоточных измерений топоцентрических дальностей до ИСЗ. Испытания выполнены для объектов с различной высотой полёта над поверхностью Земли.

Исполнение задуманного стало возможным благодаря опоре на самые основные методы классиков небесной механики и астрометрии – последовательные приближения, промежуточные орбиты, канонические преобразования:

метод последовательных приближений применялся на всех этапах, от построения моделей движения объектов в численно-аналитической форме до дифференциального улучшения параметров модели на основе наблюдений;

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

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

и, в то же время, потребовало привлечения новых идей:

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

алгоритм разложения возмущающей функции в тригонометрический ряд с численными коэффициентами является неудовлетворительным по причине существенных вариаций числовых значений эксцентриситета орбиты и угла наклонения орбиты космических тел и должен быть заменён на алгоритм преобразования возмущающего гамильтониана в сумму “элементарных” слагаемых, каждое из которых является функцией как позиционных, так и угловых параметров орбиты объекта, а также величин, характеризующих изменения возмущающих факторов модели движения;

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

Достоверность полученных результатов подтверждается примерами обработки наблюдений различных объектов, сравнением с данными, предоставляемыми Международной службой вращения Земли, и сопоставлением с параметрами движения ИСЗ, публикуемыми в Интернете.

Научная новизна работы заключается:

в совокупности формул и соотношений для вычислений в рамках обобщённой задачи двух неподвижных центров с максимально возможной методической точностью;

в методе преобразования возмущающей функции задачи, заключающейся в построении всё более сложных конструкций на основе простых начальных соотношений;

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

в совокупности алгоритмов построения моделей движения космических объектов с помощью численно-аналитического метода;

в методике расчёта ситуаций опасных сближений объектов.

Практическая ценность диссертации определяется тем, что:

предлагаемые алгоритмы численно-аналитического метода построения моделей движения справедливы в широком классе элементов орбит ИСЗ и позволяют выполнять оценку параметров движения объектов и геодинамических параметров на основе наблюдений;

алгоритмы применяются для вычисления целеуказаний на Звенигородской научной базе ИНАСАН и филиале ИНАСАН на пике Терскол; тезисы доклада “Об эфемеридном обеспечении фотометрических наблюдений ИСЗ” опубликованы в сборнике материалов пятой научно-технической конференции ФВА РВСН, с.77-78, 2000 год (совместно с Н.С.Бахтигараевым);

алгоритмы использовались для планирования космических экспериментов и предварительной редукции результатов наблюдений спутника Метеор-3М сотрудниками Центральной аэрологической обсерватории ГМЦ;

было подготовлено учебное пособие “Модель движения ИСЗ”, являющееся реконструкцией и расширением монографии профессора Е.П.Аксёнова [5].

Основные результаты, полученные в диссертации, докладывались автором и соавторами на семинарах и конференциях:

семинар по небесной механике и семинар по астрометрии ГАИШ МГУ ;

семинар ИНАСАН “Проблемы происхождения и эволюции кометноастероидного вещества в Солнечной системе и астероидная опасность” ;

конференция “Астрометрия, геодинамика и небесная механика на пороге XXI века”, Санкт-Петербург, 2000 год ;

конференция “АСТРОЭКО-2002”, Терскол, 2002 год ;

конференция “Околоземная астрономия – 2003”, Терскол, 2003 год ;

конференция “Горизонты Вселенной”, Москва, 2004 год ;

конференция “Околоземная астрономия – 2005”, Казань, 2005 год.

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

1. Герасимов И.А., Чазов В.В. Переменные действие-угол в обобщённой задаче двух неподвижных центров. //Труды ГАИШ. 1988. Т.59. С.46-52.

2. Гаипова А.Н., Чазов В.В. Комплекс программ Лента. //Измерительная техника. 1991. Т.6, С.30-30.

3. Чазов В.В. Основные алгоритмы численно-аналитической теории движения искусственных спутников Земли. // Труды ГАИШ. 2000. Т.68.

С.5-20.

4. Герасимов И.А., Чазов В.В., Рыхлова Л.В., Тагаева Д.А. Построение теории движения тел Солнечной системы, основанной на универсальном методе вычисления возмущающей функции. //Астрономический вестник.

2000. Т.34. Номер 6. С.559-566.

5. Герасимов И.А., Чазов В.В., Тагаева Д.А. Применение универсального метода вычисления возмущающей функции в численно-аналитической теории движения малых планет. //Вестник Московского университета. Серия 3. Физика. Астрономия. 2000. Номер 3. С.55-57.

6. Бахтигараев Н.С., Чазов В.В. Компьютерное моделирование условий наблюдений небесных тел. //Кинематика и физика небесных тел. 2003.

Номер 4. С.105-107.

7. Чазов В.В. Создание численно-аналитической теории движения небесных тел. /Труды конференции “Околоземная астрономия – 2003”, Терскол.

2003. С.171-175.

8. Бахтигараев Н.С., Чазов В.В. Информационное обеспечение космических экспериментов на основе численно-аналитической теории движения искусственных спутников Земли. //Космические исследования. 2005. Т.43.

Номер 5. С.386-389.

9. Бахтигараев Н.С., Чазов В.В. Моделирование движения космических аппаратов с учётом рекомендаций Международного астрономического союза. /Труды конференции “Околоземная астрономия – 2005”, Казань. 2005.

С.281-285.

10. Чазов В.В., Герасимов И.А., Соловьёва О.Д. Изотропные и гармонические координатные условия в пространстве-времени Солнечной системы.

//Вестник Московского университета. Серия 3. Физика. Астрономия.

2006. Номер 2. С.66-68.

11. Клишин А.Ф., Чазов В.В., Бахтигараев Н.С., Костюк Н.Д. Об оценке уровня техногенной опасности в зоне размещения КА “Электро-Л”. //Вопросы радиоэлектроники. Серия “Радиолокационная техника”. 2007. Выпуск 2. С.40-46.

12. Чазов В.В., Бахтигараев Н.С., Костюк Н.Д. Наблюдения спутника “Молния 3-39” в Звенигородской обсерватории ИНАСАН и определение времени падения. //Вестник СибГАУ. 2011. Выпуск 6(39). С.183-185.

13. Бахтигараев Н.С., Лёвкина П.А., Сергеев А.В., Чазов В.В. Наблюдения неизвестного фрагмента космического мусора в Терскольской обсерватории.

//Вестник СибГАУ. 2011. Выпуск 6(39). С.186-189.

На защиту выносятся следующие положения.

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

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

2. Впервые показано, что возмущающие функции различного происхождения зависят от пяти “начальных” соотношений между координатами спутника. “Начальные” соотношения с точностью, ограниченной только возможностями компьютера, были представлены в виде сумм слагаемых, зависящих от параметров промежуточной орбиты. Найден рекуррентный способ конструирования общей возмущающей функции на основе “начальных” соотношений. В предыдущих исследованиях других авторов разложения в тригонометрические ряды были получены отдельно для каждого из возмущений, обусловленных геопотенциалом, притяжением Луны и Солнца, приливными явлениями и давлением солнечного излучения.

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

4. С помощью тензорного преобразования получены приближённые формулы связи между координатами небесных тел Солнечной системы, полученными при использовании “изотропных” координатных условий с одной стороны, и “гармонических” координатных условий с другой стороны. Общий вывод состоит в следующей рекомендации: в прикладных задачах достаточно записать релятивистские уравнения движения искусственного спутника Земли на основе “гармонических” координатных условий, а при вычислении возмущающих сил использовать современные численные эфемериды Солнца, Луны и планет, полученные в метрике с “изотропными” координатными условиями.

5. Выражения для учёта возмущений от приливов упругой Земли и океанических приливов, записанные во вращающейся системе отсчёта и рекомендованные Международной службой вращения Земли, преобразованы к системе отсчёта, связанной с истинным экватором даты.

6. В результате анализа массива высокоточных лазерных измерений топоцентрических дальностей до ИСЗ Лагеос и Лагеос-2 получены оценки параметров вращения Земли и скоростей смещения измерительных пунктов в Земной системе отсчёта. Для некоторых объектов на основе наблюдений определены оценки отношения средней площади поверхности к массе спутника.

7. Разработана методика предсказания ситуаций сближения объектов искусственного происхождения в околоземном пространстве. Методика была применена на этапе проектирования параметров орбиты метеорологического “стационарного” спутника.

В первом разделе диссертации выполнена постановка задачи.

Материал второго раздела составляют основные алгоритмы модели движения ИСЗ в численно-аналитической форме.

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

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

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

Вместо процедуры разложения возмущающей функции предложен алгоритм преобразования возмущающего гамильтониана в сумму “элементарных” слагаемых. Отличительной чертой алгоритма является способ конструирования возмущающей функции на основе “начальных” функций координат.

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

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

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

Выражаю признательность профессору Нестерову Вилену Валентиновичу за простую и чёткую постановку задачи в далёком 1972 году и за предоставление возможности её решения.

Работа выполнена при финансовой и моральной поддержке гранта РФФИ № 00-02-17558, научно-образовательного проекта Института астрономии РАН и госконтракта №02.740.11.0249 “Оптический мониторинг ближнего и дальнего космического пространства роботизированной сетью телескопов МАСТЕР ” ГАИШ МГУ.

1 Постановка задачи 1.1 Стандартные соглашения Первый сборник стандарта вычислений был подготовлен в 1983 году Морской обсерваторией США в виде специального циркуляра [198]. В 1988 году, после создания Международной службы вращения Земли, был напечатан сборник “IERS Standards”. Затем последовали выпуски 1992, 1996 [188] и годов [189]. Все выпуски построены по принципу преемственности и в основных направлениях повторяют самый первый циркуляр.

Представлены факторы, которые необходимо учитывать при прогнозировании движения космических объектов. Сюда входят гравитационное поле Земли, притяжение Луны, Солнца, приливы упругой Земли и океанические приливы, давление солнечного излучения. Рекомендуется использовать современную численную модель движения Луны, Солнца и планет [202].

Были предложены новые понятия: “промежуточный экватор”, “невращающаяся начальная точка” и “угол вращения Земли” (приложение, с.161).

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

1.1.1 Системы отсчёта Международная земная опорная система отсчта (International Terrestrial Reference System, ITRS) в выбранный начальный момент времени t0 задана положениями R(t0 ) и скоростями V станций, участвующих в геодинамических исследованиях. Положение наблюдательной станции в момент времени t вычисляется по формуле где параметры смещения Ri (t) определяются на основе знаний по физике приливов. Поддержка и расширение Международной земной опорной системы координат (International Terrestrial Reference Frame, ITRF) постоянная забота Международной службы вращения Земли (International Earth Rotation Service, IERS) и Международной службы лазерной локации (International Laser Range Service, ILRS).

Международная небесная опорная система отсчта (International Celestial Reference System, ICRS) задана высокоточными экваториальными положениями 608 внегалактических радиоисточников. Точка начала координат помещена в барицентр Солнечной системы. В оптическом диапазоне эту систему представляют каталог Гиппаркос и его расширения, например, Опорный каталог Тихо [91].

В рекомендациях МАС подчёркивается, что пространство-время не может быть представлено математически в рамках единой координатной системы отсчёта. Для того, чтобы обрабатывать современные астрономические наблюдения, необходимо использовать несколько релятивистских систем отсчёта [181]. Теория преобразования координат и координатного времени, рекомендуемая МАС, представлена в работе С.М.Копейкина [86].

Небесная опорная система отсчёта с началом координат в центре масс Земли называется геоцентрической небесной системой отсчёта.

Современные численные модели движения искусственных спутников Земли построены в геоцентрической небесной системе отсчёта [100]. Модели предназначены для обработки высокоточных светолокационных наблюдений ИСЗ [103] и определения параметров вращения Земли.

Система отсчёта с началом координат в центре масс Солнечной системы получила название барицентрической небесной системы отсчёта.

В численных моделях движения Солнца, Луны и больших планет Солнечной системы, разработанных в Лаборатории реактивного движения США [202] и в Институте прикладной астрономии Российской академии наук [119], использована барицентрическая небесная опорная система отсчёта.

Промежуточное положение между двумя подходами занимает современная модель прецессии-нутации МАС2000А [183], необходимая для установления связи земной и небесной опорных систем отсчёта.

1.1.2 Пространство-время Пространство-время определяется значениями 10 компонент метрического тензора g, каждый компонент является функцией координатного времени и трёх пространственных координат, индексы и принимают значения 0, 1, 2, 3. Квадрат интервала имеет вид по повторяющимся индексам выполняется суммирование, c – скорость света, t – координатное время.

Метрические коэффициенты g находят в результате решения уравнений поля [34]. Задача имеет малый параметр: отношение скорости пробной частицы к скорости света. Отношение потенциала взаимодействия пробных частиц к квадрату скорости света пропорционально второй степени малого параметра. Решение уравнений поля в постньютоновском приближении со всей необходимой точностью соответствует достаточно медленным движениям небесных тел Солнечной системы и относительно небольшим силам их взаимного притяжения.

Симметричный метрический тензор g имеет различный вид в геоцентрической и барицентрической системах отсчёта. Десять компонент тензора связаны между собой четырьмя произвольными соотношениями, так называемыми координатными условиями.

Квадрат интервала (2) перепишем в виде Индексы i, j принимают значения 1, 2, 3. По повторяющимся индексам выполняется суммирование. Запись gij dxi dxj, например, равносильна выражению Символ Кронекера ij = 1 для i = j, и ij = 0 для i = j.

1.1.3 Окрестности Земли Компоненты метрического тензора в геоцентрической небесной системе отсчёта запишем в виде где через UE обозначено выражение для геопотенциала, а символом R обозначена сумма возмущающих функций, обусловленных действием Луны, Солнца, планет и приливов.

Известный эффект Лензе-Тирринга обусловлен действием векторного потенциала Vi. Пусть f – гравитационная постоянная, C – наибольший момент инерции Земли, z – угловая скорость вращения Земли, тогда где x, y, z – координаты точки в окрестности Земли, вычисляемые в геоцентрической системе отсчёта, Параметр t в формуле (3), получил название геоцентрическое координатное время (Geocentric Coordinate Time, TCG).

В резолюциях МАС вводится понятие земное время (Terrestrial Time, TT).

По определению, шкала TT отличается от шкалы TCG постоянным дрейфом:

где LG = 6.969290134 · 1010 является определяющей постоянной. Сама же шкала земного времени совпадает со шкалой атомного времени TAI. Различие заключается в положении нуль-пункта:

Шкала атомного времени формируется на основе показаний совокупности эталонов частоты служб времени. Всемирное координированное время UTC только на целое число секунд отличается от атомного времени TAI.

1.1.4 Солнечная система Современная теория движения планет, Луны и Солнца [202], построенная в барицентрической системе отсчта, получена численным интегрированием релятивистских уравнений движения, записанных в постньютоновском приближении с помощью изотропных координатных условий [200].

Международный астрономический союз, принимая во внимание факт, что многие работы по теории относительности выполнены при использовании “гармонических” координат, оказавшихся полезными для приложений, рекомендует выбор гармонических координатных условий [199].

В разделе 2.6.1 (с.109) представлен алгоритм учёта отличий использования изотропных и гармонических координатных условий.

Параметр t в формуле (3), получил название барицентрическое координатное время (Barycentric Coordinate Time, TCB).

Модель движения Солнца, Луны и планет, созданная в Лаборатории реактивного движения [202], скрывает в себе ещё одну особенность: аргументом для численного интегрирования уравнений движения является специально подобранная переменная Teph, шкала которой очень близка шкале земного времени TT. Эти две шкалы никогда не расходятся между собой более, чем на 2 миллисекунды. Во многих приложениях таким отличием можно пренебречь. В противоположность этому, разность между TT и TCB возрастает на 0.5 секунды в год [138].

В точке пространства, находящейся в окрестности Земли и имеющей барицентрический вектор r, разность двух шкал координатного времени TCB и TCG равна [199] где rE и vE – вектор положения и вектор скорости Земли относительно барицентра Солнечной системы, Teph – аргумент современных численных моделей движения планет, согласующая переменная LB = 1.550505 · 108.

1.1.5 Прецессия и нутация Матрицы поворота против часовой стрелки на угол вокруг осей OX, OY и OZ обозначим, соответственно, через R1 (), R2 () и R3 ().

Стандартная эпоха получила условное обозначение J2000, она совпадает с юлианской датой 2451545.0 (полдень, 1 января 2000 года).

Наклон эклиптики к экватору в стандартную эпоху J2000 обозначают символом 0. Численное значение параметра входит в список основных астрономических постоянных.

Матрица прецессии выполняет преобразование от фиксированного экватора, заданного небесной стандартной системой отсчёта, к системе подвижного экватора даты. Формулы для вычисления величин A, zA, A, получивших название прецессионных параметров Ньюкома-Андуайе, имеют вид полиномов по времени.

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

Матрица нутации образована тремя последовательными поворотами и выполняет преобразование от подвижного экватора даты к истинному где A – наклон мгновенной эклиптики к подвижному экватору.

t – время в юлианских столетиях от эпохи J2000.0.

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

Формулы для вычисления фундаментальных аргументов представляют из себя полиномы по времени.

Количество слагаемых в каждом из рядов для нутации в долготе и нутации в наклоне в модели МАС 1980 года равно 106 [188]. В модели МАС2000А ряды нутации содержат 678 лунно-солнечных и 687 планетных членов [189].

1.1.6 Стандартное преобразование Стандартное преобразование между земной (T) и геоцентрической небесной (C) системами отсчёта выполняется по формуле где P – матрица прецессии, N – матрица нутации, S – гринвичское истинное звздное время, xp, yp – координаты полюса.

Верхний штрих означает транспонирование соответствующей матрицы.

Величины являются продуктом наблюдательной деятельности. Они входят в набор параметров вращения Земли.

Гринвичское среднее звздное время определено как функция всемирного времени где TU – время в шкале всемирного времени UT1, отсчитываемое в юлианских столетиях от стандартной эпохи J2000.

Разность UT = UT1UTC всемирного времени и всемирного координированного времени входит в формулу для вычисления промежутка звздного времени от 0h UT1 до момента наблюдений, отсчитываемого от начала суток в шкале UTC:

где r = 1.002737909350795 + 5.9006 · 1011 TU 5.9 · 1015 TU.

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

1.1.7 Координаты обсерваторий В публикации Международной службы вращения Земли [189] предложен следующий алгоритм определения смещений наземных пунктов.

Координаты станции наблюдений, Луны, Солнца и сумма векторов смещения Ri вычисляются в земной опорной системе отсчёта.

Пусть f m – геоцентрическая гравитационная постоянная, f m2 – гравитационный параметр Луны, f m3 – гравитационный параметр Солнца, r,, – сферические координаты наземного пункта, r2, 2, 2 – сферические координаты Луны, r3, 3, 3 – сферические координаты Солнца, r – геоцентрический единичный вектор наземного пункта, r2 – геоцентрический единичный вектор Луны, r3 – геоцентрический единичный вектор Солнца, e – единичный вектор, направленный к востоку, n = [ e ] – единичный вектор, направленный на север:

В рамках современной теории упругости Земли, разработанной геофизиками, определены значения чисел Лява k2, k3, h2, h3 и Шида l2, l3 для второй и третьей сферических гармоник разложения приливного потенциала.

Особенность теории состоит в том, что параметры k, h, l имеют действительную и мнимую части, зависящие от частоты приливной волны. Числа h и l2 зависят от широты наземных пунктов:

Первое слагаемое суммы вычисляют по формуле Второе слагаемое суммы обусловлено третьей сферической гармоникой разложения приливного потенциала Луны:

где Зависимость от долготы учтена с помощью двух параметров:

Вклады в смещение пункта, пропорциональные l21, l22, равны где и присоединённые функции Лежандра Для мнимых частей чисел h и l предложены значения:

Вклады в смещение пункта, пропорциональные этим параметрам, равны где Зависимость действительных и мнимых частей чисел Лява и Шида от частоты приливной волны вносит свой вклад в смещение наземных пунктов:

Вариации потенциала центробежных сил вызывают следующие изменения координат обсерваторий в миллиметрах:

разность t исчисляется в годах от эпохи J2000.

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

Способ получения канонических уравнений на основе элементов промежуточной орбиты впервые дан в работе Е.П.Аксёнова [5]. В той же статье с точностью до первого порядка малости выведены формулы связи между каноническими элементами и произвольными постоянными интегрирования.

Модификацию канонических уравнений в статье [175] выполнил С.Н.Яшкин.

В монографии [5] профессор Е.П.Аксёнов отмечает: “По-видимому, канонические уравнения являются наиболее удобными для аналитических исследований. Эти уравнения были использованы С.Н.Вашковьяк для построения теории движения спутников Марса [42], Л.П.Насоновой для вычисления вековых возмущений третьего порядка в движении спутника [107] и Н.А.Сорокиным при определении долгопериодических неравенств второго порядка [129].” Описание особых, резонансных движений небесных тел является важной задачей небесной механики. С запуском искуcственных спутников Земли число таких задач пополнилось проблемой “критического” наклона. Решению их в случае близких спутников посвящены статьи С.Н.Яшкина [173] и [174].

Среди первых теорий движения “стационарных” спутников следует отметить работы М.А.Вашковьяка [36] и С.Г.Журавлва [73, 74].

С.Г.Журавлв опубликовал в монографии [76] универсальный метод исе следования острорезонансных задач небесной механики.

Особого внимания заслуживают работы А.С.Сочилиной [133, 134, 135] по эволюции орбитальных элементов геостационарных объектов и спутников на орбитах с критическими наклонениями.

Оригинальный метод учта лунно-солнечных возмущений, реализованный М.А.Вашковьяком [36], существенно улучшил А.А.Кантер [80, 81].

В цикле работ Е.Н.Поляховой подробно рассмотрены возмущения в движении спутника, обусловленные световым давлением. Прекрасный обзор результатов опубликован в специальном издании [121].

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

Оценка влияния несферичности Луны на движение ИСЗ выполнена в статьях С.Н.Вашковьяк [46] и Н.А.Сорокина [132].

Искусственные спутники Земли с высотой полёта менее 1500 километров испытывают дополнительные ускорения, обусловленные сопротивлением атмосферы и, как правило, продолжительным по времени прохождением объекта в тени Земли.

В книге Т.В.Бордовицыной и В.А.Авдюшева [24] изучены самые современные численные и аналитические методы, даны примеры решения конкретных задач и сравнительные оценки точности.

Численное интегрирование уравнений движения искусственных спутников Земли – отдельная, интересная и трудная задача. Здесь выделяются замечательные исследования М.С.Ярова-Ярового [170, 171, 172].

Заслуживает внимания оригинальный метод [130], разработанный Н.А.Сорокиным. Уравнения движения представлены в виде, позволившем объединить метод Энке и промежуточную орбиту, построенную на основе решения обобщнной задачи двух неподвижных центров.

Коллектив авторов в составе С.К.Татевян, Н.А.Сорокина и С.Ф.Зале ткина в статье [141] теоретически и практически развивает неявный одношаговый метод численного интегрирования, берущий начало от работ Е.Эверхарта [186] и Ю.В.Плахова [120].

Метод численного интегрирования уравнений движения позволяет выполнить независимую проверку результатов отдельных этапов построения аналитической теории движения небесных объектов. Об этом аспекте исследоваПостановка задачи ний смотрите, например, работы Е.П.Аксёнова и Л.М.Доможиловой [10, 11], Н.В.Емельянова [69, 70], С.М.Кудрявцева [87, 88].

Первую группу исходных данных составляют рекомендации Международного астрономического союза об использовании числовых значений астрономических постоянных и о системах отсчёта пространства-времени [199] и рекомендации Международной службы вращения Земли по составу ускорений, возмущающих движение спутника [189]. Учёт упомянутых рекомендаций необходим для того, чтобы сделать результаты обработки различных экспериментов сопоставимыми.

Для описания движения искусственных спутников приняты во внимание возмущения, обусловленные гравитационным полем Земли, притяжением Луны и Солнца, приливными деформациями упругой Земли и океаническими приливами. Учтены также влияние светового давления и эффект захода объектов в тень Земли. При создании алгоритмов использованы числовые значения коэффициентов модели гравитационного поля Земли с условным названием JGM3 [205]. Алгоритмы включают два свободных параметра, зависящих от формы и ориентации космических аппаратов – эмпирический коэффициент отражения и коэффициент эмпирического ускорения.

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

Начальные параметры движения большого числа объектов доступны пользователям в виде средних элементов кеплеровской орбиты [206]. Формат данных получил название “двустрочные элементы”. Детальное описание формата содержится в приложении на с.163.

Числовые значения топоцентрических дальностей до искусственных спутников Земли, измеренные с помощью метода лазерной локации, хранятся в виде наборов строк [190]. Каждая порция строк соответствует одному прохождению спутника в поле зрения пункта наблюдений и содержит два-три десятка “нормальных” точек, являющихся результатом процедуры осреднения большого числа “сырых” наблюдений. Точность одной “нормальной” точки находится в интервале от 2 до 10 сантиметров. Текстовый формат данных получил название “Quick Look”. Детальное описание формата содержится в приложении на с.164.

9306102 Стелла 800 0.0007 98.476 14.27121 -2.9342 0. 7501001 Старлет 900 0.0206 49.817 13.82205 3.3041 -3. 8606101 Эйджисаи 1500 0.0011 50.010 12.44446 2.5461 -3. 7603901 Лагеос 6000 0.0044 109.835 6.38665 -0.2137 0. 9207002 Лагеос-2 6000 0.0137 52.650 6.47295 0.4370 -0. 8900103 Эталон-1 19500 0.0012 65.312 2.13156 -0.0062 -0.0320 8903903 Эталон-2 19500 0.0011 64.778 2.13205 -0.0036 -0.0332 В табл.1 для нескольких спутников представлены средние значения минимальной высоты полёта hmin в километрах, эксцентриситета орбиты e, угла наклонения орбиты i в градусах, среднего движения n в оборотах за сутки, скорости изменения аргумента перигея в градусах за сутки, скорости изменения долготы восходящего узла в градусах за сутки и количество “нормальных мест” Np, обычно получаемых на интервалах наблюдений протяжённостью один месяц.

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

Методы построения модели движения космического аппарата составляют третью группу исходных данных.

Во всех центрах анализа применяют численное интегрирование дифференциальных уравнений движения [100, 204]. Численные теории имеют дело с мгновенными параметрами движения. Проблема восстановления числовых значений таких переменных на основе средних элементов орбиты не может быть решена без дополнительной информации.

Авторы работ [69, 88] представили оригинальные аналитические теории, но не сообщают о систематических испытаниях предлагаемых алгоритмов на основе наблюдений.

В данном исследовании разработаны алгоритмы и предложена модель движения ИСЗ в численно-аналитической форме [159].

1.2.1 Эволюция элементов орбиты С помощью наборов данных [206], содержащих “двустрочные элементы”, можно проследить эволюцию параметров орбиты различных объектов.

На рис.1 представлены изменения эксцентриситета и угла наклонения орбиты спутника Стелла на интервале времени 14 лет.

Эволюция элементов орбит других “геодезических” спутников из табл. происходит похожим образом. Разность максимальных и минимальных значений эксцентриситета орбиты не превосходит 0.005. Для значений угла наклонения разность всегда менее одного градуса.

1.2.2 Уравнения движения Дифференциальное уравнение в частных производных, представляющее движение частицы с массой mp в искривлённом пространстве-времени с метрическим тензором g (формулы (4)-(6) на с.20), имеет вид [93] где функция S – действие для частицы, g – контравариантный метрический тензор.

В постньютоновском приближении [34] где UE – геопотенциал, R – сумма возмущающих функций, обусловленных действием Луны, Солнца, планет и приливов.

Подстановки S S mp c 2 t и S mp S превращают выражение (17) в уравнение Гамильтона-Якоби с потенциалом UE + R и возмущающим гамильтонианом обусловленным эффектами общей теории относительности.

Вектор V определён на с.20. Вектор скорости спутника в геоцентрической небесной системе отсчёта обозначен символом v.

Пусть в геоцентрической небесной системе отсчёта задан вектор r с координатами x, y, z. Параметр t соответствует шкале всемирного координированного времени. Точка над буквой означает дифференцирование по параметру t. Тогда Дифференциальному уравнению в частных производных (18) эквивалентна система трёх дифференциальных уравнений второго порядка:

где параметр C – наибольший момент инерции Земли, z = 0.7292115 · радиан в секунду – угловая скорость вращения Земли.

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

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

Выбранная неинерциальная система вращается относительно небесной опорной системы отсчёта с угловой скоростью (t). С точностью до первых степеней малых величин компоненты вектора (t) равны При этих ограничениях уравнение Гамильтона-Якоби (18) эквивалентно системе канонических уравнений движения [93] где r – вектор в подвижной системе отсчёта, p = r+ r – обобщённый импульс, K – гамильтониан задачи, 1.2.3 Геопотенциал Пусть f m – геоцентрическая гравитационная постоянная;

r0 – экваториальный радиус Земли;

Jn – числовые коэффициенты при зональных гармониках;

Cnk, Snk – числовые коэффициенты при тессеральных и секториальных гармониках разложения гравитационного поля Земли в ряд по сферическим функциям.

В земной опорной системе отсчёта выражение для геопотенциала вне поверхности Земли имеет вид где r,, – геоцентрические расстояние, широта и долгота;

Pn (sin ) – полиномы Лежандра;

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

В статье большого коллектива авторов [205] опубликованы числовые значения постоянных f m, r0, Nmax, Cnk, Snk, соответствующих модели JGM3.

Нормированные коэффициенты Cnk, Snk и числовые коэффициенты формулы (24) связаны соотношениями:

Коэффициент C20 по порядку величины близок значению 1.0·103. Такой же порядок величины имеют как коэффициент при второй зональной гармонике J2, так и параметр, характеризующий сжатие (сплюснутость) земного эллипсоида. Численные значения нормированных коэффициентов Cnk и Snk при n = 2 находятся на уровне 1.0 · 106. Принято говорить, что коэффициент J2 является величиной первого порядка малости относительно сжатия. В этой терминологии все остальные коэффициенты разложения гравитационного поля Земли имеют второй порядок малости относительно сжатия, а релятивистское слагаемое Kr является величиной четвёртого порядка малости.

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

Через r и v обозначим вектор положения и вектор скорости космического аппарата относительно центра Земли в системе экватора и эклиптики, фиксированных на стандартную эпоху J2000.0 (геоцентрическая небесная опорная система отсчёта, с. 17).

FE ускорение, обусловленное геопотенциалом, FM ускорение, вызываемое притяжением Луны, FS ускорение, вызываемое притяжением Солнца, Ft ускорение, обусловленное приливами упругой Земли, Fo ускорение, обусловленное океаническими приливами, Fp ускорение, вызываемое притяжением больших планет, Fr ускорение, обусловленное давлением солнечного света, Fa ускорение, обусловленное торможением в атмосфере.

Уравнения движения имеют вид Ускорение в движении искусственного спутника, вызываемое внешним телом (Луной, Солнцем или планетой) с массой mp, имеет вид где r – геоцентрический вектор объекта, r = |r| – модуль геоцентрического вектора объекта, rp – геоцентрический вектор планеты, rp = |rp | – модуль геоцентрического вектора планеты.

Первое слагаемое – это главная часть ускорения, второе слагаемое, не зависящее от положения спутника, называют косвенной частью [165].

Формуле (26) соответствует потенциал также составленный из двух частей, главной и косвенной. Учитывая неравенство r < rp и разлагая в ряд величину, обратную расстоянию между спутником и планетой, получим формулу где Pn (z) – полином Лежандра порядка n. Суммирование начинается со ника и не принимается во внимание. Слагаемое ряда при n = 1 взаимно уничтожается с косвенной частью.

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

Алгоритм вычисления векторов FE, Ft и Fo в небесной системе отсчёта включает в себя дифференцирование выражения для геопотенциала UE (с.35), заданного в земной системе отcчёта. Последовательность формул, составляющих этот алгоритм, представлена в приложении на с.167.

Для определения ускорения, обусловленного давлением солнечного излучения, вычисляют вектор положения спутника r(t) rS (t ) относительно Солнца с учтом времени распространения света. Ускорение имеет вид где aE – астрономическая единица, c – скорость света, Cr – коэффициент отражения (с.81).

Сопротивление верхней атмосферы Земли в движении спутников учитывается с помощью аэродинамического ускорения, направленного противоположно вектору относительной скорости объекта [64]:

где Sb = CD – баллистический коэффициент, CD – аэродинамический коэффициент сопротивления, A, m – средняя площадь поверхности и масса космического аппарата, (h) – плотность атмосферы на высоте h, единица измерения – кг/м3, e – вектор угловой скорости вращения Земли, h = r r0 1 2 – высота полёта спутника в километрах над поr верхностью общеземного эллипсоида, r = |r| – геоцентрическое расстояние объекта, z – координата объекта по оси аппликат, r0 = 6378.1363 км – экваториальный радиус Земли, = – сжатие Земли, безразмерный параметр.

Пусть средняя площадь поверхности аппарата задана в квадратных метрах, а масса – в килограммах, тогда для согласования единиц измерений априорное значение величины Sb следует вычислять по формуле Плотность атмосферы является сложной функцией параметров, характеризующих солнечную активность и геомагнитную обстановку в атмосфере Земли. Для вычисления значения плотности верхней атмосферы Земли рекомендована модель ГОСТ Р 25645.166-2004 [17].

1.2.5 Приливный потенциал Первый член суммы (27) является основным в теории приливов:

Пусть rp, p, p – сферические координаты Луны или Солнца, вычисляемые в земной опорной системе отсчёта, тогда Теорема сложения для полиномов Лежандра позволяет записать где Функции разлагают в ряды. Член ряда определён численным значением амплитуды и аргументом тригонометрической функции, синус или косинус. Аргумент заm висит от линейной комбинации среднего звёздного времени S и фундаментальных аргументов lM, lS, FM, D, M. Аргумент имеет фазу и частоту.

Приливы вызывают дополнительный потенциал где k20, k21, k22 – числа Лява. В современной теории приливов эти параметры состоят из действительной и мнимой частей и зависят от частоты приливной волны, совпадающей с частотой, которую имеет аргумент соответствующего слагаемого рядов (32).

Подставляя (31), (32) в (33) и сравнивая с (24), получим где k20, k21, k22 – номинальные значения приливных чисел Лява, не зависящие от частоты приливной волны.

Аналогичную процедуру применяют и для второго члена суммы (27), сферической гармоники третьего порядка:

где k30, k31, k32, k33 – числа Лява.

Рекомендовано вычислять малые вариации числовых значений коэффициентов C40, S40, C41, S41, C42, S42 по формуле Зависимость действительных и мнимых частей параметров k20, k21, k от частоты приливной волны рекомендуется учитывать с помощью формул [189]:

где Численные значения амплитуд и коэффициентов содержатся в специальных таблицах [189].

Динамические эффекты, обусловленные океаническими приливами, предлагается учитывать таким же образом:

где Cnk, Snk, nk – амплитуды и аргумент приливной волны порядка n и степени k, имеющей частоту [185], Fnk – численные коэффициенты.

Пусть f = 6.673 · 1011 м3 /(кг · с2 ) – гравитационная постоянная, ge = 9.7803278 м/с2 – ускорение свободного падения на экваторе, w = 1025 кг/м3 – плотность морской воды, тогда где kn – коэффициенты нагрузочной деформации, Потенциал центробежных сил изменяется по причине движения полюсов.

Деформации Земли, вызванные этим эффектом, приводят к появлению дополнительного потенциала. Влияние его предложено учитывать с помощью поправок [189] где текущие значения координат полюса xp, yp измеряются в секундах дуги, разность t исчисляется в годах от эпохи J2000.

1.3 “Стационарный” объект “Стационарные” спутники совершают приблизительно один оборот за одни земные сутки. Многие из них являются управляемыми объектами. С помощью специальных двигателей несколько раз в год выполняется коррекция среднего движения и угла наклонения орбиты спутника. В результате таких действий объект продолжительное время находится в экваториальной плоскости Земли и с точностью до нескольких угловых минут имеет заданную географическую долготу. Когда управление параметрами движения объекта прекращают, начинается эволюция элементов орбиты под действием гравитационных сил и сил, вызываемых давлением солнечного излучения.

Одним из таких объектов является COMSTAR 4 (номер 81018A по каталогу КОСПАР). Он был запущен на “стационарную” орбиту в 1981 году и лет отработал в “точке стояния” над восточным побережьем Южной Америки. После окончания активной стадии работы объект был переведён в область орбиты с географической долготой около 75.

На основе средних элементов орбиты в формате “NORAD” (с.163) вычислим начальный вектор состояния спутника. На интервале времени 55 лет выполним численное интегрирование дифференциальных уравнений движения с учётом ускорений, обусловленных геопотенциалом, притяжением Луны и Солнца и световым давлением.

На рис. 2 (с.44) представлены графики изменения оскулирующих элементов орбиты и географической долготы спутника, вычисленных в момент прохождения восходящего узла. Один из периодов изменения эксцентриситета орбиты близок к 11 годам, а минимальное и максимальное числовые значения эксцентриситета различаются более чем в 5 раз. В случае с углом наклонения ситуация ещё более интересная: в продолжение 27 лет этот элемент орбиты убывал от 15 градусов до почти нулевого значения, а затем начал возрастать. Интересен также следующий эффект: сдвиг по фазе в вариациях большой полуоси и географической долготы.

Эволюция угла наклонения орбиты показывает большие отличия как по амплитуде, так и по периодам колебаний в сравнении с соответствующими значениями для “геодезических” спутников (рис.1 на с.32).

Амплитуды колебаний большой полуоси и эксцентриситета орбиты “стационарных” объектов сравнимы с соответствующими значениями для “геодезических” спутников (рис.1 на с.32). Эволюция угла наклонения орбиты показывает большие отличия как по амплитуде, так и по периодам колебаний. Для описания движения искусственных спутников Земли, находящихся на “стационарных” орбитах, необходима разработка специальных моделей.

1.4 Промежуточный потенциал В аналитических теориях движения ИСЗ выражение (24) (с.35) представляют в виде суммы двух функций:

промежуточного потенциала W и возмущающей функции RE.

42170. ``` ` ``` ``` ``` ```` `` `` ``` `` `` ``` ` ``` ` `` `` ```` ```` ``` ``` ```` `````` ` ``` ````` ``` ```` ` ` ``` `` ` `````` ```` ` `` ``` `` ` ```` ````` `` ```` ` ` `` `` ```` ` ``` ` `` ` `` `` ``` `` ``` ```` ````` `` ` `` ` ``` ` ``` ` ````` `` ```` `` `````` ``` ` ` ` ` ` ` ` ````` ``` `` ```` ` `` `` ``` ```` ``` ``` `` ```` `` ` ` ``` ` ` `` ```` ``` ``` `` ` ``` ` ```` ```` ``` `` ``` `` ` `` `` ` `` ````` ``` ```` `` `` ```` ````` ` ` ``` ``` ``` ``` ``` ``` `` ` ``` ``` ` ` `` ` ` ``` ``` `` ```` ```` ```` ```` ``` ```` ```` `` ``` ` `` `` `` ` ` `` `` `` ``` ` ` `` ``` ``` ` ` ` ` `` ` `` ```` ```` `` ```` ``` `` ``` `` `` ` `` ``` ``` ```` `` ``` ``` ``` ` ``` ``` ` ` ````` `` ``` ` `` ``` `` ``` ```` ``` `` `` `` `` `` `` ` ` ``` ` ` ```` ```` ````` ```` `` ` ``` ` ` ` `` ` `` ```` ```` ``` ```` `` ` ```` ``` `` ` ` ` `` ` `` ` ` ``` `` ` ``` `` `` ` ```` ` ` ```` ``` ```` ```` `` ```` ``` ``` `` `` ` ` ` ``` ` `` ``` ` `` ```` ``` ````` ````` `` ``` ` ` ````` 42160. 0. 0. 14.720 T 0. 80.143 `` `` `` ``` `` `` `` `` `` ``` ` `` ``` `` `` ``` `` ``` `` ``` `` `` `` `` ``` ` `` `` `` `` `` ``` `` `` `` ``` `` ``` `` ``` ``` `` `` ``` `` `` `` `` `` `` ` `` `` `` ``` ```` ```` ```` ``` ``` ``` ``` ``` ```` ``` ``` ``` ``` ``` `` ``` ``` ``` ```` ```` ``` ``` ``` ``` ``` ``` ``` ``` ```` ```` ``` ``` ```` ```` ``` ```` ```` ``` ``` `` ``` ``` ```` ``` ``` ``` ``` ```` ``` ``` ``` ``` ```` 75. 69. Потенциал W выбирают так, чтобы, с одной стороны, он включал в себя самые существенные по величине слагаемые, и, с другой стороны, уравнение можно было проинтегрировать. Общий интеграл уравнения и шесть произвольных постоянных определяют промежуточную орбиту. Неравенства, обусловленные влиянием как возмущающей функции RE, так и других факторов, вычисляются с помощью теории возмущений [2].

Если положить и использовать сферические координаты, то уравнение Гамильтона-Якоби (39) можно проинтегрировать. Решением является кеплеровская промежуточная орбита.

Кеплеровская промежуточная орбита послужила основой для фундаментальных исследований А.М.Фоминова [154, 155]. Численно-аналитический метод для вычисления возмущений кеплеровской орбиты впервые применил М.А.Вашковьяк [36, 37]. Этот же подход для решения своей задачи использовала И.В.Тупикова [147, 148].

где c и – вещественные постоянные, Если положить Промежуточный потенциал (41) – (43) впервые был предложен Е.П.Аксёновым, Е.А.Гребениковым и В.Г.Дёминым [7]. Потенциал W может рассматриваться как силовая функция задачи двух неподвижных центров с комплексно-сопряжёнными массами (m/2) (1+ 1) и (m/2) (1 1), удалнными друг от друга на мнимое расстояние, равное 2 1 c.

Задачу с такой силовой функцией называют обобщнной задачей двух неподвижных центров или задачей Гредеакса [98].

После перехода к сжатым сфероидальным координатам,, w, которые связаны с прямоугольными координатами x, y, z формулами обобщённая задача двух неподвижных центров интегрируема в квадратурах.

Обращение квадратур впервые выполнил Е.П.Аксёнов в работе [4].

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

Теория движения ИСЗ, использующая промежуточную орбиту, основанную на решении обобщённой задачи двух неподвижных центров, была построена Е.П.Аксёновым и изложена в монографии [5]. Результаты обработки наблюдений содержатся в статьях [8, 9].

Более детальная по количеству принимаемых во внимание возмущающих факторов аналитическая теория была разработана К.В.Холшевниковым и Е.В.Тимошковой [157, 142, 143].

Формулы для учёта влияния возмущающих факторов с точностью до третьего порядка относительно малого параметра – сжатия Земли – были получены Н.В.Емельяновым [68, 69] и В.А.Тамаровым [139].

1.5 Канонические преобразования С помощью обобщённых координат r и обобщённых импульсов p запишем канонические уравнения в неинерциальной системе отсчёта (выражение (22) на с.35) в виде:

В векторных уравнениях (46) K0 является гамильтонианом обобщённой задачи двух неподвижных центров, а возмущающая функция K2 состоит из слагаемых, имеющих второй порядок малости относительно сжатия (формула (23)):

Интегрирование уравнений с гамильтонианом K0 изложено в монографии профессора Е.П.Аксёнова [5]. Система имеет интеграл энергии. Произвольная постоянная интегрирования, связанная с интегралом энергии, была обозначена символом 1. Характеристическая функция K0 была приравнена произвольной постоянной 1. После перехода к сжатым сфероидальным координатам (45) было составлено уравнение Гамильтона – Якоби:

Полный интеграл уравнения был найден методом разделения переменных:

где 3 – произвольная постоянная, 1 и 1 – постоянные, а многочлены четвёртой степени () и F () имеют вид где 2 – ещё одна, третья произвольная постоянная.

Общий интеграл будет даваться уравнениями где 1, 2 и 3 – произвольные постоянные, и С помощью новой независимой переменной, связанной с временем t дифференциальным соотношением и общего интеграла задачи (52), (53) можно записать две эллиптические квадратуры где c3 и c4 – постоянные интегрирования, и два эллиптических интеграла, первый для вычисления долготы а второй для установления неявной зависимости переменных t и В уравнениях (56) и (57) c5 и c6 – произвольные постоянные, а t0 – начальный момент времени.

Далее в монографии [5] канонические уравнения (46) в переменных r и p были преобразованы в канонические уравнения с использованием канонических элементов L, G, H, l, g, h промежуточной орбиты, построенной на основе решения обобщённой задачи двух неподвижных центров:

где K = 1 + K2 – гамильтониан системы.

В случае отсутствия возмущений K2 = 0, и система (58) имеет решение причём скорости изменения угловых переменных являются постоянными величинами.

Отметим следующий факт: необходимо знать только частные производные от гамильтониана K по каноническим элементам, представлять же гамильтониан в виде явной функции переменных действия L, G, H нет никакой необходимости.

Формулы для вычисления канонических элементов L, G, H получены в монографии [5] с точностью до первого порядка относительно сжатия. Модификацию канонических уравнений в статье [175] выполнил С.Н.Яшкин.

Приближённое решение уравнений (58) может быть получено с помощью метода канонических преобразований [156].

Выполним переход от оскулирующих элементов L, G, H, l, g, h к новым переменным L, G, H, l, g, h. Функцию преобразования и вековой гамильтониан обозначим, соответственно, S2 и K2. Эти функции, как и возмущающий гамильтониан K2, имеют второй порядок малости относительно сжатия.

В новом гамильтониане K2 ( t, L, G, H, l, g, h ) оставим те слагаемые возмущающей функции K2, период изменения которых превышает несколько суток. Функцию преобразования S2 найдём по формуле С точностью до второго порядка относительно малого параметра, сжатия Земли, выражения S2, K2, K0 = 1, K2 можно считать функциями как оскулирующих элементов L, G, H, l, g, h, так и новых переменных L, G, H, l, g, h.

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

1.6 “Осреднённые” уравнения Уравнения движения, записанные в новых переменных, назовём “осреднёнными” уравнениями. Они сохраняют каноническую форму и имеют вид:

1.7 Постановка задачи Применение численно-аналитического метода построения модели движения искусственных спутников Земли заключается:

• в разработке алгоритма вычисления параметров промежуточной орбиты, основанной на решении обобщённой задачи двух неподвижных центров, с точностью, ограниченной только возможностями компьютера;

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

• в разработке алгоритма преобразования возмущающей функции в сумму слагаемых, зависящих от параметров промежуточной орбиты, с точностью, ограниченной только возможностями компьютера;

• в разработке алгоритмов аналитического интегрирования слагаемых возмущающей функции и выделения слагаемых “осреднённого” гамильтониана;

• в разработке алгоритмов численного интегрирования “осреднённых” уравнений движения и учёта короткопериодических неравенств;

• в определении параметров модели на основе наблюдений.

2 Основные алгоритмы 2.1 Промежуточная орбита В монографии Е.П.Аксёнова [5] получены формулы в буквенном виде.

Зависимость между переменными в этих формулах – явная. Точность вычисления промежуточной орбиты ограничена вторым порядком малости относительно сжатия Земли. Исследования Е.П.Аксёнова, Н.В.Емельянова и В.А.Тамарова [12] показали, что формулы удовлетворяют многим практическим приложениям. С помощью специальной программы “универсальный пуассоновский процессор” [31], созданной коллективом под руководством В.А.Брумберга, Н.В.Емельянов [67] довл решение до четвртого порядка оте е носительно малого параметра – сжатия Земли.

В данной работе предлагается другой подход [159]. С помощью неявных соотношений между переменными будет построен алгоритм для проведения расчётов с точностью, ограниченной только возможностями компьютера. Совокупность формул, выражающая как явные, так и неявные зависимости между переменными, становится алгоритмом для решения поставленной задачи и может быть запрограммирована. Входными данными вычислительной процедуры будут числовые значения величин f m, c, и начальный вектор состояния. Результатом станут численные значения параметров промежуточной орбиты.

2.1.1 От вектора состояния к параметрам орбиты Пусть для начального момента времени t = t0 известны прямоугольные координаты x = x0, y = y0, z = z0 и скорости x = x0, y = y0, z = z0.

Требуется определить элементы промежуточной орбиты a, e, и угловые переменные l0, g0, h0 на момент t0.

Сначала по известным координатам и их производным находим постоянные интегрирования 1, 2, 3. Затем определяем группу позиционных элементов a, e,. После чего найдём начальные значения l0, g0, h0.

2.1.2 Постоянные интегрирования 1, 2, Введм следующие обозначения [5] Для момента t = t0 имеем где через 0, 0, w0 обозначены сфероидальные координаты для начального момента времени.

Начальные значения производных от сфероидальных координат по времени t равны Численные значения произвольных постоянных интегрирования определены формулами [5] 2.1. Методом итераций на основе уравнений где находим числовые значения величин Обозначим = +. Для вычисления поправки используем неявное уравнение в котором Параметр 2 является величиной первого порядка малости относительно сжатия. При вычислениях по формулам (68) – (73) с точностью, ограниченной возможностями компьютера, достаточно выполнить 7 итераций.

Определим также числовые значения величин 1 и Далее последовательно найдём числовые значения параметров Величины 2, e, k2, k2 определены выражениями Вспомогательные величины n, n равны где Параметры 1, s, k1, k1, d, определены выражениями Вспомогательные величины m, m равны где Числовое значение параметра равно а параметры 2j, 2j определены выражениями где K(k) – полный эллиптический интеграл первого рода с модулем k, Cn – “биномиальные коэффициенты”:

Разложение для K(k) имеет вид Ряд сходится при k < 1.

Параметр является малой величиной первого порядка относительно сжатия. Параметры 2j, 2j пропорциональны сжатию в степени j, их числовые значения быстро убывают при увеличении индекса j. По этой причине верхний предел знака суммирования в формулах (81) и (82), равный бесконечности, можно с учтом вычислительной точности заменить на некоторое целое число M. Коэффициент при второй зональной гармонике разложения геопотенциала имеет порядок 103, и численное значение M = 7 обеспечит точность вычислений с 20 значащими цифрами.

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

2.1.4 Алгоритм операций с полиномами Дифференциальные соотношения между независимой переменной и угловыми переменными, имеют вид где параметры k1 и k2 являются величинами первого порядка малости относительно сжатия.

Для дифференциальных соотношений справедливы разложения где верхний предел суммирования, равный бесконечности, заменён на конечное число M = 7. Правые части равенств (87) и (88) можно рассматривать как многочлены от переменных cos и cos, соответственно.

Неявная связь между угловыми переменными и задаётся формулой Это выражение и формулы (80), (81) (с.56) для вычисления параметров, 2j, 2j были получены в результате приравнивания дифференциалов (87) и (88) и интегрирования с использованием формул При определении числовых значений параметров (85) и (86) в дополнении к равенствам (90) используем ещё два полезных алгоритма.

Первый алгоритм опубликован в монографии Н.В.Емельянова [66].

Пусть заданы многочлены переменной z с численными коэффициентами ai, bj. Результатом умножения полиномов будет многочлен с численными коэффициентами ck Интересен и алгоритм деления многочленов или, другими словами, алгоритм выделения целой части. Пусть I > J, тогда Вычисление коэффициентов dk и gn происходит за I J + 1 шагов. Начальные значения dk = ak, 0 k I и gn = 0. Далее на каждом шаге изменяются численные значения J коэффициентов dk, причм коэффициент dk при текущей старшей степени z k становится равным нулю, и степень многочлена в числителе уменьшается на единицу. Это достигается вычитанием из текущего многочлена в числителе многочлена в знаменателе, умноженного на выражение Параметры 0, 0, µ, j, j появляются как результат рассмотрения эллиптических интегралов [5]:

Для вычисления первого интеграла учтём разложение (87) и соотношение где Применяя алгоритм умножения полиномов (91), вычислим многочлен относительно cos, являющийся произведением правой части равенства (87) и выражения 1 2d cos + d2 cos2. Выделяя далее с помощью алгоритма (92) целую часть, получаем На основе соотношений (90) вычислим интеграл где 0, i – численные коэффициенты.

Во втором интеграле (93) выполним разложение подытегрального выраc жения в ряд по степеням Отношение с помощью равенства может быть представлено в виде многочлена по степеням cos. Действительно, где q – величина первого порядка малости относительно сжатия. Применяя к правой части равенства алгоритм умножения многочленов, получаем После возведения (97) в квадрат получим искомый многочлен с численными коэффициентами Выполняя M 1 раз операцию умножения текущего многочлена на выражение (98) и складывая, получим многочлен для подынтегральной функции (95). Умножим его на дифференциал (88) и проинтегрируем, учитывая формулы (90). В результате оказывается, что Определяя параметры придадим формуле для вычисления переменной вид Параметры, j, j появляются как результат рассмотрения эллиптических интегралов [5]:

Первый интеграл запишем в виде Параметры p и p2 + q 2 вычисляются по формулам (78) (с.55) и являются величинами первого порядка малости относительно сжатия.

Для вычисления интегралов учтм, что и представим последовательно выражения в виде полиномов по степеням cos.

Далее с помощью алгоритма умножения полиномов составим подынтегральное выражение и проинтегрируем его. Умножая результат на числовой множитель, получим Для вычисления второго интеграла заметим, что 2 можно представить как произведение полиномов а для дифференциала d использовать формулу (87).

Дважды применяя алгоритм умножения многочленов, после интегрирования получим Определяя параметры, ip, if формулами придадим окончательный вид равенству, связывающему неявным образом угловые переменные l,, и E :

2.1.5 Угловые переменные l0, g0, h Найдм значения угловых переменных E,, в начальный момент t = t0. Переменная E (аналог эксцентрической аномалии):

Переменная (аналог истинной аномалии):

Переменная (аналог аргумента широты):

Численные значения угловых переменных l, g, h в начальный момент времени t0 определяются с помощью соотношений 2.1.6 От параметров орбиты к вектору состояния Пусть в произвольный момент времени t известны численные значения позиционных элементов a, e, и угловых переменных l, g, h. Требуется вычислить координаты x(t), y(t), z(t) и скорости x(t), y(t), z(t).

По формулам (68) – (72) найдём численные значения произвольных постоянных интегрирования 21, 2, 3. Учтём, что для углов наклонений орбиты, превышающих 90, постоянная 3 < 0, то есть sign(3 ) = 1.

Далее с помощью выражений (73) – (84) и алгоритмов, представленных в разделе 2.1.4 (с.57 – с.62), определяем числовые значения параметров и коэффициентов. Методом последовательных приближений найдм значения уге ловых переменных E,, :

используя при этом соотношения (103).

Переменную вычисляем по формуле Сфероидальные координаты равны где выражение для переменной w следует из равенств Производные по времени от сфероидальных координат равны Искомые прямоугольные координаты определены формулами:

Искомые производные по времени от прямоугольных координат определены формулами:

2.1.7 Проверка алгоритмов В модели гравитационного поля Земли JGM3 геоцентрическая гравитационная постоянная f m, экваториальный радиус Земли r0, коэффициенты при второй J2 и третьей J3 зональных гармониках имеют следующие числовые значения [205]:

Подставляя в (42) (с.45) вместо r0, J2 и J3 их числовые значения (114), получим c = +209.729971234944476 км, В табл.2 содержатся численные значения нескольких коэффициентов, соответствующих модели гравитационного поля Земли JGM3 [205] и промежуточному потенциалу (формула (41) на с.45):

Таблица 2: Числовые значения коэффициентов зональных гармоник J4 = 1.6193312052 · 106 J4 = 1.166177041845644980 · J5 = 0.2277161017 · 106 J5 = +0.005469555070724556 · J6 = +0.5396484905 · 106 J6 = +0.001249751229132727 · J7 = 0.3513684419 · 106 J7 = 0.000008844878268002 · J8 = 0.2025187153 · 106 J8 = 0.000001332336307242 · J9 = 0.1193687135 · 106 J9 = +0.000000012692302667 · J10 = 0.2480568648 · 106 J10 = +0.000000001412746230 · J11 = +0.2405652138 · 106 J11 = 0.000000000017045753 · J12 = 0.1819117030 · 106 J12 = 0.000000000001489618 · J13 = 0.2075677323 · 106 J13 = +0.000000000000021939 · J14 = +0.1174173876 · 106 J14 = +0.000000000000001561 · J15 = 0.0176272697 · 106 J15 = 0.000000000000000027 · J16 = +0.0311943085 · 106 J16 = 0.000000000000000002 · Для иллюстрации утверждения о том, что точность алгоритмов вычисления промежуточной орбиты ограничена только накоплением ошибок округления, был проделан численный эксперимент.

Для модельного объекта, совершающего 13.98 оборотов за сутки и имеющего следующие параметры орбиты методом численного интегрирования дифференциальных уравнений движения на интервале времени 100 суток были вычислены эталонные значения прямоугольных координат xe (t), ye (t), ze (t). В правых частях уравнений движения были учтены гравитационные силы, обусловленные только промежуточным потенциалом (табл.2).

Начальному вектору состояния объекта на момент t = 0. x = +4917.49973747459503 км, x = 1.2636786137103486 км/с, y = +3693.31783253124247 км, y = +6.0704892431019494 км/с, z = +3866.34490247898799 км, z = 3.9703600780539020 км/с соответствуют параметры промежуточной орбиты (с.52) На моменты t эталонных значений координат с помощью набора параметров (116) и алгоритма (с.63) были определены положения модельного объекта xm (t), ym (t), zm (t). Величина биты. График изменения r(t) представлен на рис.3. Отношение величины нения, равного 8.64 · 106 секунд времени, не превышает 1012.

2.1.8 Новые результаты В разделе 2.1 представлены следующие новые результаты, полученные автором впервые [159].

• Выведена система уравнений (68)-(73), предложен и реализован итерационный алгоритм её решения.

• Предложен и реализован алгоритм вычисления эллиптических интегралов на основе операций над многочленами (87)-(102).

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

2.2 Частные производные 2.2.1 Переменные действия Переменные действия A1, A2, A3 определены формулами [5] где () и F () даются равенствами (51), 1 и 2 – корни многочлена (), между которыми изменяется координата, а 1 и 2 – корни многочлена F (), лежащие в промежутке [1, +1].

Вторая группа элементов – переменные угол B1, B2, B3, сопряжнные переменным действия A1, A2, A3, определены формулами [5] где S – выражение (50), в котором 1, 2, 3 считаются функциями A1, A2, A3.

Три уравнения общего интеграла (52) (с.48) можно записать в виде Сравнивая выражения (50), (52) и (117), сразу получаем формулы для частных производных Учитывая пределы интегрирования, имеем Обращение матрицы приводит к частным производным от параметров k по элементам Ak :

Выражения (117) перепишем в виде Подставляя сюда явный вид функций (), F () (51) и сравнивая результат подстановки с соотношениями (52), замечаем, что почти все интегралы, кроме двух новых, нам уже известны:

2.2.2 Частные производные по элементам L, G, H Канонические элементы действия L, G, H определены формулами [5] то есть где, µ, 0, 0,, I, I – величины первого порядка малости относительно сжатия, причём Численные значения I и I находятся с помощью алгоритмов умножения и интегрирования многочленов относительно cos и cos соответственно (с.57), после подстановки пределов интегрирования число в знаменателе сокращается.

Частные производные от постоянных интегрирования 1, 2, 3 по переменным действия L, G, H принимают вид Первое полезное свойство совокупности частных производных (120) проявляется в алгоритме вычисления наборов величин a, e, и 21, 2, на основе известных численных значений элементов L, G, H. Для удобства обозначим L1 = L, L2 = G, L3 = H. В нулевом приближении полагаем Далее на основе уравнений (68)-(73) методом последовательных приближений находим числовые значения величин (a)0, (e)0, ()0, а затем, также в нулевом приближении, последовательно вычисляем другие параметры задачи: элементы (L)0, (G)0 и частные производные от величин 1, 2 по этим элементам. Уточннные значения параметров находим по формуле линейных разностей Процесс сходится очень быстро, за две – три итерации.

Второе полезное свойство соотношений (120) тоже достаточно очевидно. С их помощью выполняется дифференцирование позиционных параметров a, e, по каноническим элементам действия L, G, H. Алгоритм состоит из трх шагов. Продифференцируем каждую из формул (68), (69), (70) по явно входящим величинам a, e2, 2 (приложение, с.171) и вычислим матрицу Численное обращение этой матрицы даст значения частных производных от параметров a, e2, 2 по 1, 2, 3, и, наконец, численные значения производных по каноническим элементам получаем в результате суммирования 2.2.3 Частные производные по параметрам a, e, Частные производные первого порядка по каноническим элементам действия L, G, H будут необходимы при вычислении правых частей осреднённых уравнений движения (раздел 2.5, c.102). Промежуточная орбита на основе обобщённой задачи двух неподвижных центров включает в себя несколько десятков параметров (списки под номерами (76), (85), (86), например), зависящих от позиционных элементов a, e, явным образом. Если известны производные от любого параметра p по a, e,, то производные от этого параметра по L, G, H вычисляются с помощью соотношений (122):

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

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

Основу алгоритма составляет множество M одномерных массивов чисел.

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

Для каждого параметра p промежуточной орбиты определён свой массив чисел P(p). Первый элемент массива содержит численное значение параметра, следующие элементы содержат численные значения трёх первых и шести вторых независимых частных производных по позиционным элементам:

Численные значения десяти независимых частных производных третьего порядка по параметрам a, e, составляют содержание остальных элементов массива:

Для двух массивов чисел P(p), Q(q) определены операции сложения и умножения. Пусть a1 = a, a2 = e, a3 = и индексы i, j, k принимают значения 1, 2, 3. Результатом сложения P(p) + Q(q) является массив чисел R(r), образованный попарным сложением элементов:

Результатом умножения P(p) Q(q) является массив чисел R(r) :

Для операции возведения в произвольную степень b, выполняемой над массивом чисел P(p), достаточно определить значения элементов массива R(r) :

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

Алгоритм формирования конечного множества M заключается в вычислении массивов, содержащих частные производные от величин 2 (формула (71) на с.54), Q (формула (72)), 2 1 (формула (68)), 2 (формула (69)), 3 (формула (70)), и всех остальных, включая упомянутые ранее параметры из списков под номерами (76), (85), (86). Алгоритм составлен так, что сумма или произведение параметров заменяются на операции сложения или умножения, выполняемые над массивами чисел, соответствующих этим параметрам. Деление равносильно возведению в степень, равную (1).

2.2.4 Производные от угловых переменных Частные производные от угловых переменных,, E, по позиционным параметрам a, e, получаются путём непосредственного дифференцирования формул (103)-(106):

и использования метода последовательных приближений. В нулевом приближении частные производные определены формулами:

Достаточно двух-трёх итераций.

Частные производные от переменных,, E, по каноническим угловым элементам l, g, h вычисляются по той же схеме:

с производными в нулевом приближении 2.2.5 Новые результаты В разделе 2.2 представлены следующие новые результаты, полученные автором впервые [159].

• Разработан алгоритм вычисления численных значений канонических переменных действия.

• Разработан алгоритм вычисления частных производных по каноническим элементам промежуточной орбиты.

• Разработан алгоритм вычисления частных производных высших порядков по позиционным параметрам.

• Точность вычислений ограничена только возможностями компьютера.

Присутствие частных производных высших порядков позволяет сохранить информацию о том, как тот или иной параметр был вычислен. Пусть, например, одномерный массив чисел, соответствующий параметру p, получен для конкретных числовых значений позиционных элементов a0, e0, 0. Расчёты показали, что для малых численных значений разностей мгновенное значение любого параметра p с относительной точностью порядка 104 107 определяется с помощью первых членов формулы Тейлора:

где 2.3 Алгоритм преобразования В публикациях А.М.Фоминова [154, 155], В.В.Нестерова и Г.В.Романовой [108, 111] авторы использовали различные промежуточные орбиты (кеплеровскую орбиту и орбиту, построенную на основе решения К.Акснеса [176]), но разложение возмущающей функции было выполнено с помощью формул кеплеровского движения, то есть с точностью до второго порядка относительно малого параметра – сжатия Земли.

В статьях Н.В.Емельянова [68], Н.В.Емельянова и Л.П.Насоновой [71] с помощью некеплеровской промежуточной орбиты, основанной на решении обобщённой задачи двух неподвижных центров, получено разложение возмущающей функции, обусловленной гравитационным полем Земли, притяжением Луны и Солнца, с точностью до третьей степени относительно сжатия.

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

2.3.1 Аномальный геопотенциал Пусть x, y, z – координаты вектора r в системе отсчёта, связанной с истинным экватором.

Модуль расстояния равен r = | r | = Координаты вектора r в земной системе отсчёта обозначим x, y, z.

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

Величина в этом выражении – угол долготы в земной системе отсчёта.

Угол долготы, отсчитываемый от истинной точки весеннего равноденствия вдоль истинного экватора, равен w = S +, причём S – гринвичское истинное звёздное время (с.24), и Разложим выражение для полинома Лежандра Pn (/r) в ряд Тейлора.

С точностью до первого порядка относительно z (формула (143)) получим Аномальную часть геопотенциала в системе истинного экватора запишем в виде где величины Cnk = Cnk cos kS Snk sin kS, Snk = Cnk sin kS + Snk cos kS, являются функциями звёздного времени.

Принимая во внимание определение присоединённых функций Лежандра (формула (237) на с.172 приложения), формулу (144) и две теоремы сложения для суммы аргументов тригонометрических функций (239), получим следующий результат: рекуррентный алгоритм расчёта возмущающего потенциала (145) основан на начальных значениях четырёх функций координат спутника Знаменатель, возникающий при вычислении cos kw, sin kw, сокращается с таким же выражением в числителе после раскрытия соотношеk) ния для присоединённой функции Лежандра Pn (z/r) (приложение, с.172).

2.3.2 Притяжение светил Возмущающий потенциал, обусловленный притяжением Луны или Солнца, можно представить в виде (27) (с.37) в разложении оставлено Np 1 слагаемых, индекс p в случае Луны принимает значение M, в случае Солнца – S, rp – геоцентрический вектор положения возмущающего тела, rp – геоцентрическое расстояние, Угол долготы возмущающего тела, отсчитываемый от истинной точки весеннего равноденствия вдоль истинного экватора, обозначим wp. Раскроем выражение и применим теорему сложения для полиномов Лежандра (формула (238), с.172 приложения):

Формула (27) выглядит теперь так где коэффициенты зависят только от координат возмущающего объекта Сравнение выражений (145) и (147) показывает, что рекуррентный алгоритм расчёта потенциала (147), обусловленного притяжением третьего тела, основан на начальных значениях четырёх функций координат спутника 2.3.3 Эффект приливов Гравитационное действие Луны и Солнца вызывает приливы упругой Земли. Во внешнем пространстве появляется дополнительный потенциал (формула (33) на с.40) где knk – числа Лява. В современной теории приливов параметры k2k состоят из действительной kR и мнимой kI частей и зависят от частоты приливk 2k ной волны. Параметры k3k являются действительными числами. Коэффиp) (p) циенты Ank, Bnk зависят от координат возмущающего объекта в системе истинного экватора и определены соотношениями (148).

Рекуррентный алгоритм расчёта дополнительного потенциала (150) основан, как и в случае геопотенциала (145), на начальных значениях четырёх функций координат спутника (146) (с.79) :

2.3.4 Световое давление Выражение для ускорения Fr, обусловленного давлением солнечного света, записано на с.38 (формула (28)).

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