Содержание к диссертации
Введение
1 Методы теории специальных возмущений в задачах околопланетной динамики 16
1.1 Линеаризация и регуляризация 17
1.1.1 Основные принципы методов линеаризации и регуляризации 17
1.1.2 Системы уравнений Шперлинга-Боде и Кустаанхеймо-Штифеля
1.2 Сглаживающие преобразования 19
1.3 Численная стабилизация
1.3.1 Неустойчивость кеплеровского движения 21
1.3.2 Диссипативный метод Баумгарта 23
1.3.3 Консервативный метод Баумгарта 26
1.3.4 Стабилизация по времени 27
1.3.5 Метод Накози 27
1.3.6 Стабилизация в случае почти кругового движения
1.4 Метод вариации постоянных 28
1.5 Метод Энке
1.5.1 Основные принципы метода Энке 30
1.5.2 Классический метод Энке 30
1.5.3 Уравнения Энке в KS-переменных 31
1.5.4 Метод Энке для приведения систем к стандартному виду 32
1.5.5 Улучшение опорной орбиты
1.6 Проблема короткопериодпческих возмущений 36
1.7 Сравнительный анализ эффективности методов
1.7.1 Численный эксперимент 40
1.7.2 Выбранные объекты 41
1.7.3 Интегратор Эверхарта 42
1.7.4 Характеристики эффективности численного интегрирования 42
1.7.5 Численные результаты 43
1.8 Интегратор Гаусса-Эверхарта для численного решения
дифференциальных уравнений первого порядка 49
1.8.1 Основные формулы 50
1.8.2 Интегрирование на шаге 52
1.8.3 Формулы интегратора как одно из представлений неявного метода Рунге-Кутты 53
1.8.4 Повышение порядка интегратора 53
1.8.5 Выбор шага з
Методы сглаживания орбит с короткопериодическими возмущениями в задачах динамики далеких спутников планет 56
2.1 Возмущенный гармонический осциллятор 58
2.2 Ограниченная круговая задача трех тел 59
2.3 Преобразования уравнений 60
2.4 Усредненные уравнения 61
2.5 Уравнение энергии и его решение в круговом случае 62
2.6 Модифицированные усредненные уравнения движения 64
2.7 Приближенная численная оценка остаточных вековых эффектов 66
2.8 Моделирование притяжения гауссова кольца
в пространственном случае 67
2.9 Численные результаты 69
Гравицентрические системы координат в задачах спутниковой и астероидной динамики 73
3.1 Ошибки округления в кеплеровских членах. Задача двух тел 74
3.2 Ограниченная задача трех тел. Астероидная задача 76
3.3 Временные ошибки 78
3.4 Метод синхронного слежения 81
3.5 Исследование эффективности интегрирования при использовании гравицентрических координатных систем 83
3.6 Влияние методических ошибок в положении планеты на точность определения движения сближающегося астероида 86
Методы решения обратных задач динамики близких спутников. особенности и их исследование 89
4.1 Проблема неоднозначного определения орбит 90
4.2 Круговая задача 91
4.3 Обратная задача орбитальной динамики 100
4.4 Методы решения обратной задачи
4.4.1 Метод Гаусса-Ньютона 101
4.4.2 Демпфированный метод Гаусса-Ньютона 102
4.4.3 Метод Левенберга-Марквардта 102
4.4.4 Метод Гельфанда-Цетлина 103
4.4.5 Составной метод 103
4.5 Задача двух тел. Исследование эффективности методов для решения обратной задачи 104
4.6 Модель спутникового движения 111
4.7 Наблюдения спутников 115
4.8 Определение орбитальных параметров 117
4.9 Другие оценки 122
4.10 Сравнение с эфемеридами JUP230 126
4.11 Условие неоднозначности в определении орбит близких спутников 128
5 Методы оценивания точности орбитальных параметров в обратных задачах околопланетной динамики. особенности и их исследование 130
5.1 Методы статистического моделирования возможных значений параметров для исследования точности орбит 130
5.1.1 Задача наименьших квадратов и доверительные области 131
5.1.2 Метод возмущенных наблюдений 132
5.1.3 Метод возмущенных оценок 134
5.1.4 Случай неравноточных наблюдений 137
5.1.5 Бутстрэп метод 138
5.1.6 Численный пример
5.2 Общее описание ошибок в определяемых спутниковых орбитах 145
5.3 Построение областей возможных значений параметров для новых спутников Юпитера 148
5.3.1 Размеры начальных вероятностных областей 148
5.3.2 Вероятностные области через оборот 149
5.3.3 Зависимость между временными интервалами наблюдаемости спутников и размерами вероятностных областей 152
5.3.4 Нелинейное оценивание 154
5.3.5 S/2003 J02: спутник или астероид? 157
5.4 Особенности в оценивании точности орбитальных параметров для близких спутников 160
5.5 Быстрые отображения для статистического оценивания вероятности попадания объекта в малый объем 164
5.5.1 Несингулярный случай 165
5.5.2 Сингулярный случай 167
Заключение 170
Список литературы
- Сглаживающие преобразования
- Ограниченная круговая задача трех тел
- Ограниченная задача трех тел. Астероидная задача
- Демпфированный метод Гаусса-Ньютона
Введение к работе
Появление новых астрометрических средств наблюдения за последние десятилетия вызвало небывалое повышение точности и стремительное увеличение количества наблюдательной информации о движении как уже известных, так и постоянно открываемых небесных тел, что к настоящему моменту естественным образом ставит перед специалистами теоретической астрономии актуальную проблему о пересмотре существующих и разработке новых математических моделей, интерпретирующих наблюдательный материал. В связи с широким использованием в астрономии вычислительных технологий математические модели сейчас все чаще реализуются на компьютерах посредством соответствующих численных методов. Таким образом, для обработки современных наблюдений от численно реализуемых математических моделей требуется беспрецедентно высокая эффективность.
Эффективность численного моделирования характеризуется точностью вычисляемых результатов, а также быстродействием их получения. Центральными результатами моделирования орбитального движения являются векторы динамического состояния исследуемого небесного тела: в прямых задачах это — прогнозируемые векторы на заданные моменты времени, получаемые на основе начальных, тогда как в обратных это — начальные векторы динамического состояния, получаемые на основе сопоставления результатов моделирования прямой задачи с наблюдениями. В этом смысле построение численной модели по наблюдениям можно рассматривать как процесс совместного решения прямой и обратной задач, что предусматривает использование самых разнообразных методов для повышения точности и быстродействия вычислений. Говоря об эффективных методах численного моделирования, обычно имеют в виду такие методы, которые повышают либо точность при сохранении быстродействия на несколько порядков и выше, либо быстродействие при сохранении точности в несколько раз и более, либо, в лучшем случае, обе характеристики эффективности одновременно.
Эффективность численного моделирования определяется следующими основными факторами: 1) решаемой задачей (особенностью значений параметров, их точностью, налагаемыми условиями); 2) принятой математической моделью для ее решения; 3) методом компьютерной реализации модели; а также 4) архитектурой процессора, на котором выполняются вычисления, и ошибками компьютерной арифметики. Трудности моделирования, вызываемые как раз первыми двумя из перечисленных факторов, ставят сейчас перед небесными механиками ряд насущных проблем, для разрешения которых необходимы конструктивные подходы.
Уникальность параметрических значений в исследуемой задаче определяет различные особенности моделируемой динамики: конфигурацию орбиты, регулярность и скорость движения небесного тела и т.д., что, разумеется, непосредственно влияет на эффективность моделирования, основанного на численном интегрировании дифференциальных уравнений как формализации орбитальной динамики. В этом случае эффективность численного моделирования может быть повышена путем выбора более удачной формализации орбиты. Различные варианты дифференциальных уравнений движения в небесной механики предоставляют так называемые методы теории специальных возмущений, преобразующие классические уравнения к лучшему виду с точки зрения их численного интегрирования, что в результате позволяет повысить методическую и вычислительную точность, а также быстродействие численного моделирования.
Период развития методов теории специальных возмущений насчитывает уже несколько столетий. Первые методы были предложены еще в XIX веке (например, метод Энке или метод вариации произвольных постоянных Лагранжа). По понятным причинам изначально они разрабатывались и применялись в теориях общих возмущений. Бурное развитие и появление оригинальных методов теории специальных возмущений приходится на 60-70 гг. прошлого века, что, очевидно, было вызвано широким использованием в небесно-механической практике компьютерных технологий. Несмотря на то, что с тех пор прошло довольно много времени и по методам вышло множество работ, обзоры которых можно найти в (Херрик, 1977; Рой, 1981; Aarseth, 1988), фактически не было ни одной, где бы проводился сравнительный анализ разнотипных методов и давались четкие рекомендации по их использованию.
Численный подход к решению задач орбитальной динамики остро ставит проблему интегрирования короткопериодических возмущений. Эта проблема связана с таким явлением, когда в правых частях дифференциальных уравнений движения присутствуют быстроосциллирующие члены, частота которых существенно превышает среднее движение небесного тела. Даже если эти члены по величине меньше кеплеровского, задаваемого центральной силой, быстродействие численного интегрирования, тем не менее, будет определяться именно короткопериодическими членами, несмотря на то, что их вклад в решение может быть ничтожно малым. Как нетрудно заметить, подобные особенности возникают при численном интегрировании плохообусловленных (жестких) дифференциальных уравнений (Miranker, 1981).
Проблема короткопериодических возмущений в небесной механике типично имеет место в задачах, где моделируется орбита медленного (внешнего) объекта под гравитационным влиянием быстрых (внутренних), как, например, в задачах о движении далеких астероидов и спутников (Авдю-шев, 2006а; Авдюшев, 2007), где учитываются влияния соответственно от планет земной группы и от близких массивных спутников типа галилеевых.
В случае малого вклада короткопериодических возмущений в исследуемую орбиту для разрешения вышеуказанной проблемы, как правило, прибегают к огрублению формализации движения путем включения возмущающих масс в центральную. Если же вклад короткопериодических возмущающих сил достаточно большой с точки зрения представления наблюдательных данных, обращаются к усреднению Гаусса (Дубошин, 1961; Jacobson, 2000; Emelyanov, 2005), что также приводит к упрощению динамической модели, но она оказывается все же точнее, нежели ее аналог с увеличенной центральной массой. Впрочем, второй подход имеет серьезную брешь, выраженную в том, что усреднение Гаусса не учитывает остаточные вековые эффекты в орбитальной долготе, которые на длительных интервалах времени могут привести к большим отклонениям моделируемых положений небесного тела относительно действительных.
В численном исследовании астероидных орбит наблюдается существенная потеря вычислительной точности в положении объекта при каждом его тесном сближении с большой планетой. Для разрешения этой трудности обычно прибегают к различного рода регуляризирующим преобразовани-
ям дифференциальных уравнений движения (Stiefel, Scheifele, 1971; Heggie, 1974; Aarseth, Zare, 1974). Однако недостаток такого подхода главным образом связан со значительным усложнением формализации орбитального движения, в то время как существуют другие более простые способы повышения вычислительной точности (Авдюшев, 2000; Авдюшев, 2003b), которые изложены в содержательной части данной диссертационной работы.
Параметры орбитальной модели, как правило, определяются из наблюдений в рамках задачи наименьших квадратов. Если ошибки наблюдательных данных носят случайный характер, то параметрическая точность будет зависеть не только от точности наблюдений, но также и от их количества: при одинаковой дисперсии случайных ошибок наблюдений точность параметров повышается с увеличением объема измерительной информации. Важным этапом в определении орбитальных параметров является оценивание параметрической точности, что обычно выполняется в контексте линейной задачи наименьших квадратов. В то же время использование линейных оценок не всегда обосновано, в особенности когда ошибки параметров достаточно большие, поскольку связь между орбитальными параметрами и представляемыми моделью положениями небесного тела, вообще говоря, нелинейна.
Обычно используемый на практике стандартный подход для оценивания параметрической точности состоит в построении границы доверительной области, определяющей вероятностный разброс параметрических ошибок, в соответствии с заданной изоповерхностью целевой функции (Press et al., 1987). Однако, если даже часто изоповерхность достаточно хорошо описывает границу доверительной области, то указать вероятность накрытия ею реальных значений параметров в нелинейном случае оказывается невозможным (Draper, Smith, 1981).
Для нелинейного оценивания параметрической точности вообще нет универсальных методов (Draper, Smith, 1981; Демиденко, 1981; Bates, Watts, 1988), однако специфика исследуемой задачи порой все же позволяет сконструировать подходящие методы, пригодные для практического использования. Подобные методы в небесной механике предлагаются, например, в работах (Milani, 1999; Virtanen et al., 2001; Muinonen et al., 2006). Эффективный подход для нелинейного оценивания реализован (впрочем, без всякого обоснования) на сайте , который вполне при-
емлем в тех случаях, когда внутренняя нелинейность (Beale, 1960) подпространства оцениваемых параметров в пространстве измеряемых равна нулю, хотя в практических задачах это условие не выполняется.
Параметрические ошибки могут быть вызваны не только ошибками наблюдений, но и особенностями обратной задачи, имеющей атрибуты некорректности. Обратная задача орбитальной динамики, как правило, сводится к минимизации некоторой целевой функции, которая выражает степень близости наблюдаемых и моделируемых положений объекта. Вообще говоря, совершенно нет никаких оснований полагать, что целевая функция имеет единственный минимум. Напротив, в большинстве обратных задач благодаря главным образом периодичности орбитального движения гиперповерхность целевой функции в пространстве определяемых параметров имеет довольно сложную структуру с многочисленными минимумами, что очень часто имеет место в нелинейных обратных задачах. Обычно не возникает затруднений в выборе минимума, соответствующего наилучшим оценкам орбитальных параметров, и таким минимумом является абсолютный, в котором значение целевой функции существенно меньше, нежели в других минимумах. Однако ситуация с выбором не всегда складывается столь благоприятно: в некоторых задачах (Авдюшев, Баныцикова, 2008) минимумы целевой функции могут быть почти равнозначными и тогда критерий качества по абсолютному минимуму не может рассматриваться как безусловный. Проблема множественности минимумов часто имеет место в обратных задачах близких спутников и, насколько нам известно, она еще здесь совершенно не изучена.
Кроме того, в обратных задачах орбитальной динамики близких спутников существует также проблема поиска решения. Дело в том, что ввиду сложности целевой функции используемые для поиска решения квазиньютоновские итерационные методы типа Гаусса-Ньютона порой оказываются неэффективными: последовательные приближения искомых параметров, получаемые итерационно, часто медленно сходятся, либо расходятся вовсе, и тогда следует прибегать к более изощренным подходам.
Цели и задачи исследования
С учетом сформулированных выше проблем при выполнении данной диссертационной работы была поставлена цель разработать методы для повышения эффективности численного моделирования орбит, исследовать возможности методов в задачах околопланетной орбитальной динамики, а также получить с их помощью результаты, имеющие прикладное значение для динамической астрономии.
Для достижения поставленной цели были решены следующие задачи.
-
Разработаны оригинальные методы теории специальных возмущений для повышения методической и вычислительной точности, а также быстродействия численного моделирования (интегрирования) орбит. В частности, предложен метод сглаживания для эффективного интегрирования астероидных орбит, а также методы так называемой стабилизации по времени и стабилизации почти кругового движения. Кроме того, получены различные эффективные модификации метода Энке применительно к интегрированию орбит в пространстве Кустаанхеймо-Штифеля.
-
Теоретически обоснованы предложенные и уже существующие и широко используемые на практике методы теории специальных возмущений применительно к характерным задачам орбитальной околопланетной динамики.
-
Проведен численный эксперимент с целью исследовать возможности методов теории специальных возмущений в задачах динамики спутников и астероидов, имеющих тесные сближения с большими планетами; выполнен сравнительный анализ эффективности методов; выработаны рекомендации по их использованию.
-
Исследована проблема короткопериодических возмущений в рассматриваемых задачах и разработаны эффективные подходы для ее решения на основе огрубления численных моделей.
-
Теоретически обосновано применение так называемых гравицентри-ческих координатных систем для повышения вычислительной точности численного интегрирования в астероидных задачах, исследованы
особенности их применения, а также получены оценки эффективности применения соответствующих алгоритмов.
-
Исследованы особенности обратных задач динамики близких спутников Юпитера, получены оценки их орбитальных параметров по имеющимся наблюдениям объектов, а также изучена состоятельность оценок с точки зрения их использования для адекватного моделирования спутниковых орбит.
-
Разработаны эффективные алгоритмы для решения обратных задач динамики близких спутников на основе комплексного использования известных итерационных методов теории оптимизации, а именно Гаусса-Ньютона и градиентного спуска, совместно с так называемым проекционным методом.
-
Разработан эффективный метод для нелинейного оценивания параметрической точности динамических моделей. Метод основан на имитации статистики Фишера и реализуется путем многократного решения нелинейной задачи наименьших квадратов при различных выборках генерируемых наблюдений, получаемых путем внесения соответствующих случайных возмущений в моделируемые положения небесного тела. Оценены параметрические неопределенности в моделируемых орбитах близких и новых далеких спутников Юпитера.
Объект и предмет исследования
Объектом исследования в работе является орбитальная околопланетная динамика, тогда как предметом исследования — методы, предназначенные для повышения эффективности численного моделирования околопланетной динамики. Таким образом, в ходе выполнения диссертационной работы предполагалось разработать новые и усовершенствовать уже существующие методы для повышения эффективности численного моделирования околопланетной динамики, исследовать их возможности, а также получить с их помощью результаты, имеющие прикладное значение для динамической астрономии.
Наш интерес к задачам околопланетной динамики вызван, прежде всего, следующими обстоятельствами: 1) спутниковые системы планет предоставляют широкий спектр орбит, на примере которых можно всесторонне исследовать возможности рассматриваемых методик; 2) движение в околопланетном пространстве с точки зрения численного моделирования усложнено рядом характерных особенностей, а именно высокими скоростями близких спутников и астероидов, тесно сближающихся с планетой; гравитационным влиянием Солнца на движение далеких спутников, приводящим к значительным изменениям их орбит, а также короткопериодиче-скими гравитационными возмущениями от близких массивных спутников типа галилеевых; 3) за последнее время в спутниковых системах открыто много новых объектов, что естественным образом привлекает к ним особое внимание. Наконец, отметим, что несмотря на принятые ограничения, применяемые методики могут быть также весьма полезны и для численного решения иных задач, не рассматриваемых в работе, например, астероидных или планетных, которые в плане моделирования имеют тесное родство с задачами околопланетной динамики.
Научная новизна исследования
Новизна диссертационной работы может быть охарактеризована следующими результатами, полученными соискателем.
-
Предложены новые методы и оригинальные модификации существующих методов теории специальных возмущений для повышения эффективности численного моделирования околопланетных орбит.
-
Проведен численный эксперимент по выявлению возможностей представленных методов применительно к задачам околопланетной динамики и на основе полученных численных результатов выполнен сравнительный анализ эффективности методов, а также выработаны четкие рекомендации по их использованию.
-
Сформулирована и исследована проблема короткопериодических возмущений, возникающая при численном интегрировании орбит далеких (медленных) объектов под гравитационном влиянием близких
(быстрых); предложены оригинальные способы эффективного решения проблемы.
-
Выявлена одна из важнейших причин существенной потери вычислительной точности при интегрировании орбит астероидов на моменты их тесных сближений с большими планетами; предложен так называемый метод синхронного слежения для повышения точности моделирования астероидного движения во время тесных сближений.
-
Сформулирована и детально исследована проблема неоднозначного определения орбит в обратных задачах динамики близких спутников.
-
Представлен оригинальный составной подход для эффективного численного поиска решения обратной задачи в случае моделирования орбит близких спутников.
-
Построена новая высокоточная численная модель спутникового движения применительно к спутниковой системе Юпитера на основе новой версии интегратора Гаусса-Эверхарта в редакции автора диссертации.
-
Используя построенную модель, получены новые оценки орбитальных параметров для всех новых далеких и внутренних спутников Юпитера.
-
Предложен новый подход типа Монте-Карло для нелинейного оценивания неопределенностей в орбитальных параметрах, получаемых из наблюдений.
10. Впервые получены оценки точности орбитальных параметров для новых спутников Юпитера.
Следует заметить, что результаты по пп. 5, 7, 8, 10 были получены при совместном участии Баныциковой М.А., где однако ее вклад имеет отношение лишь к экспериментальной части работы, которая в то же время выполнялась под руководством автора данной диссертации. Впрочем, нельзя не отметить, что по п. 5 соавтором экспериментально была обнаружена
важная с прикладной точки зрения особенность в обратных задачах динамики близких спутников, проявляющаяся в неоднозначном определении спутниковых орбит по нескольким группам наблюдений, тогда как впоследствии автор диссертационной работы дал теоретическое объяснение этому явлению.
Результаты, представленные в диссертации, также опубликованы в ряде статей с другими соавторами, однако их вклад выпадает за рамки научной работы соискателя.
Практическая значимость работы
Представленные в работе методы, а также разработанное на их основе программно-математическое обеспечение может быть использовано для повышения эффективности численного моделирования орбит в околопланетном пространстве, например, с целью идентификации и планирования наблюдений небесных тел. Как уже отмечалось, применяемые методики вполне приемлемы и для численного исследования иных задач, не рассматриваемых в работе, которые в плане моделирования имеют сходство с задачами околопланетной динамики. В частности, методы нелинейного оценивания параметрической точности могут быть весьма полезными в задачах астероидной опасности для оценки вероятности столкновения объектов с Землей, особенно в тех случаях, когда астероидная орбита определяется по немногочисленным наблюдениям на очень короткой дуге и потому имеет большие параметрические ошибки.
Некоторые предложенные автором методики, имеющие отношение к прямым задачам орбитальной динамики, применялись для решения ряда прикладных задач динамической астрономии в рамках следующих работ:
-
НИР по базовому финансированию «Исследование движения, распределения и эволюции орбит малых тел Солнечной системы по наблюдениям с Земли и из космоса» (ТГУ 3.5.96Ф, 1996-2000 гг., рук. проф. Т.В. Бордовицына);
-
НИР по базовому финансированию «Математическое моделирование движения, распределения и орбитальной эволюции малых тел
Солнечной системы по результатам измерений» (ТГУ 3.4.01Ф, 2001 2005 гг., рук. проф. Т.В. Бордовицына);
-
НИР по базовому финансированию «Исследование динамики больших популяций малых тел Солнечной системы» (ТГУ 1.36.06Ф, 2006-2010 гг., рук. проф. Т.В. Бордовицына);
-
Федеральная целевая научно-техническая программа Минпромнауки РФ «Исследования и разработки по приоритетным направлениям развития науки и техники» по теме «Исследование миграции малых тел Солнечной системы и развитие методов обнаружения потенциально опасных небесных тел, включая фрагменты космического мусора» (40.022.1.1.1108, 2002-2004 гг.);
-
Федеральная целевая научно-техническая программа Минпромнауки РФ «Исследования и разработки по приоритетным направлениям развития науки и техники», раздел «Фундаментальные исследования в области физических наук» по теме «Исследование возможностей использования радиоастрономической сети «КВАЗАР-КВО» для решения задач астрогеодинамики и фундаментального координатно-временного обеспечения России» (40.022.1.2.1109, 2002-2004 гг.).
Апробация результатов исследования
По результатам исследований опубликовано более 40 работ, из которых 5 в зарубежных изданиях (Bordovitsyna et al., 1997; Avdyushev, Bordovitsyna, 2000; Titarenko et al., 2000; Bordovitsyna et al., 2001; Avdyushev, 2003c); 12 в российских изданиях (Авдюшев, 1999а; Авдюшев, 2003а; Авдюшев, 2003b; Авдюшев, 2004; Авдюшев, 2006а; Авдюшев, 2006b; Авдюшев, Бордовицына, 2006; Баныцикова, Авдюшев, 2006b; Авдюшев, 2007; Авдюшев, Баныцикова, 2007b; Авдюшев, Баныцикова, 2008; Авдюшев, 2009), входящих в перечень ведущих рецензируемых научных изданий, в которых должны быть опубликованы основные научные результаты диссертации на соискание ученой степени доктора наук; а также 15 статей в других изданиях (Авдюшев, 1997; Авдюшев, 1998; Бордовицына и др., 1998а; Бордовицына и др., 1998b; Авдюшев, 1999b; Васильченко, Авдюшев, 1999; Авдюшев,
2000; Бордовицына, Авдюшев, 2001; Баныцикова, Авдюшев, 2002; Мишкин и др., 2002; Козаногина, Авдюшев, 2002; Авдюшев, 2005; Авдюшев, 2006с; Баныцикова, Авдюшев, 2006а; Авдюшев, 2008).
Во всех работах, выполненных, в частности, совместно с Т.В. Бордо-вицыной, вклад соавторов касался либо постановочной части, либо экспериментальной, что, впрочем, не имеет непосредственного отношения к теме данной диссертации. В других совместных работах соавторы принимали участие лишь в проведении численных экспериментов, некоторые результаты которых (а именно полученные совместно с М.А. Баныцико-вой) приводятся в диссертации. Между тем вклад соискателя как соавтора непосредственно касался разработки методов моделирования орбитального движения, причем к этой части работы другие соисполнители не привлекались.
Результаты исследований также докладывались и обсуждались на 22 конференциях различного уровня. Кроме того, некоторые результаты прошли апробацию экспертов по работам, финансируемым РФФИ (96-02-17999-а; 98-02-16491-а; 01-02-17266-а; 02-02-06888-мас; 05-02-17043-а; 08-02-00359-а) и КЦФЕ (ЕО2-11.0-6), в которых соискатель участвовал и как соисполнитель, и как руководитель.
Основные положения и результаты, выносимые на защиту
Соискатель выносит на защиту следующие положения и результаты.
-
Новые методы и оригинальные модификации существующих методов теории специальных возмущений для повышения эффективности численного моделирования околопланетных орбит.
-
Результаты сравнительного анализа эффективности методов теории специальных возмущений в задачах околопланетной динамики, а также рекомендации по их использованию.
-
Предложенные соискателем способы для эффективного решения проблемы короткопериодических возмущений.
-
Обоснование к применению гравицентрических координатных систем в астероидных и спутниковых задачах.
-
Результаты исследования проблемы неоднозначного определения орбит в обратных задачах динамики близких спутников.
-
Оригинальный составной подход для эффективного численного поиска решения обратной задачи в случае моделирования орбит близких спутников.
-
Новый подход типа Монте-Карло для нелинейного оценивания неопределенностей в орбитальных параметрах, получаемых из наблюдений.
Сглаживающие преобразования
Для приведения в соответствие координат х с физическим временем t систему урав нений необходимо пополнить уравнением Следует заметить, что для повышения точности численного интегрирования можно, как и выше, заменить уравнение времен ного преобразования на уравнение временного элемента, хотя тогда систему необходимо будет очередной раз расширить и дополнить уравнением энергии h. ( В небесно-механической практике широко используются такие сглаживающие преобразования, где в качестве / выбираются следующие величины: х (эксцентрическая аномалия); х3/2 (эллиптическая аномалия); х2 (истинная аномалия); х-1 (дуга орбиты) (Brumberg, 1992).
В астероидных задачах, где исследование орбиты астероида усложняется заметным влиянием от больших планет, например, при тесных сближениях, можно использовать следующее преобразование (Авдюшев, 2006b):
Как уже отмечалось, ляпуновская неустойчивость дифференциальных уравнений при численном интегрировании создает благоприятные условия для культивирования всевозможных ошибок, неизбежно сопровождающих любой численный процесс. Ошибки на текущем шаге интегрирования становятся ошибками начальных данных следующего, которые в дальнейшем усиливаются неустойчивостью шаг за шагом. Поэтому устойчивые уравнения для численного интегрирования более предпочтительны.
Задача стабилизации заключается как раз в том, чтобы ослабить (или в лучшем случае устранить) влияние ляпуновской неустойчивости на численное решение и улучшить таким образом поведение неустранимых ошибок интегрирования.
Родоначальниками стабилизации в небесной механике бесспорно можно считать И. Ба-умгарта (Baumgarte, 1972; Baumgarte, 1973) и П. Накози (Nacozy, 1971). Их стабилизирующие методы основаны на применении известных интегралов, которые содержат дополнительную информацию о решении и рассматриваются как ограничения, налагаемые на решение. Исторически так сложилось, что методы стабилизации не были удостоины должным вниманием со стороны небесных механиков, несмотря на то, что их исключительные возможности открываются именно в задачах небесной механики и прежде всего в задачах исследования долговременной эволюции орбитального движения. В результате ниша эволюционных задач была занята симплектическими интеграторами, которые по сути обладают достоинствами стабилизирующих методов, но при этом значительно сложнее для практической реализации. Вместе с тем, как только появились первые стабилизирующие методы, они стали широко применяться в классической механике, где и продолжают развиваться до настоящего времени.
Тем не менее, конечно же, нельзя не отметить некоторые не связанные друг с другом работы после Баумгарта и Накози (Murison, 1989; Avdyushev, 2003с; Fukushima, 2003; Zhu, Wu, 2007), в которых авторы раскрывают достоинства численной стабилизации применительно к задачам моделирования орбитальной динамики.
Технически стабилизация достигается путем исправления численного (ошибочного) решения за его уклонение от интегральной поверхности в фазовом пространстве интегрируемых переменных. В возмущенной задаче, когда интегральная поверхность динамична, а ее интегральный параметр становится переменным, для оценки уклонений решения используются те же интегральные соотношения, однако система уравнений дополняется уравнением для интегрального параметра, которое интегрируется численно совместно со всей системой (Baumgarte, 1973; Бордовицына, Сухоплюева, 1980).
Стабилизация применяется как непосредственно к самому решению в процессе интегрирования (Nacozy, 1971), так и посредством введения дополнительных так называемых стабилизирующих членов в дифференциальные уравнения (Baumgarte, 1973). Очевидно, в этом случае стабилизация тем эффективнее, чем слабее возмущения и чем медленнее меняются параметры опорных интегральных поверхностей.
Из всех интегральных соотношений авторы стабилизирующих методов выделяют энергетические, поскольку, как показывает практика, именно стабилизация по энергии является наилучшим подспорьем в борьбе с ляпуновской неустойчивостью.
Пользуясь формулами задачи двух тел, нетрудно показать (Avdyushev, 2003с), что если средние движения двух близких кеплеровских решениіі отличаются на величину An в некоторый начальный момент времени t0, при t — to — со это приводит к расхождению решений Дх и Дх в соответствии с приближенной оценкой
Оценка (1.14) показывает, что эволюция Дх носит периодически вековой характер: расхождение положений на орбитах со временем растет по величине, то локально увеличиваясь к перицентру, когда скорость возрастает, то локально уменьшаясь к апоцентру, когда скорость убывает. Также видно, что средняя скорость роста Ах явно зависит от относительного отклонения в среднем движении: чем меньше Ап/п, тем ниже скорость Дх. Однако какой бы малой ни была An, всегда наступит такой момент времени, когда Дх станет недопустимо большой. Фактически это говорит о том, что кеплеровское движение неустойчиво по Ляпунову. Впрочем, если для двух орбит An = 0, то расхождение положений Дх будет ограничено. Этот факт как раз принимается во внимание для построения эффективных методов стабилизации, применяемых в численном интегрировании орбит.
Если одно из двух решений рассматривать как точное, а другое — ошибочное, обусловленное начальными ошибками Дхо и Дхо, которые вызывают ошибку в среднем движении An, то формулу (1.14) можно использовать как оценку ошибок Дх и Дх на любой другой момент времени. При численном интегрировании кеплеровской орбиты малые ошибки типа Дхо и Дхо возникают на каждом шаге интегрирования, вследствие чего ошибка в среднем движении растет почти линейным образом (рисунок 1.1).
Ограниченная круговая задача трех тел
Интересно заметить, что величина вековых эффектов в быстрой переменной главным образом зависит от отношения размеров орбит а. второго и третьего тел и их начальной конфигурации, задаваемой углом сро. Причем для любого а 1 можно подобрать такие іро, при которых вековые эффекты будут отсутствовать, т.е. 5 = 0. Так, при достаточно малых а корни (р0 будут близки к значениям (р0 = 7г/4 + ітг/2 (г = 0,..., 3), которые фактически являются корнями уравнения cos 2 /?о = 0.
Мы исследовали точность представления движения уравнениями (2.22) в той же планетной задаче, что мы рассмотрели выше. На рисунке 2.3 показаны расхождения решений уравнений (2.22) и (2.9) в зависимости от величины а. Здесь график этой зависимости обозначен как (1 + 5)(fJ,s + [ р + А4)- Интегрирование выполнялось на интервале времени 84 оборота третьего тела10. Для сравнения мы также приводим оценки расхождений решений для невозмущенных уравнений: С(а) — 0 (fis) , для уравнений, где масса второго включена в массу первого: С (а) = 1, Со(а) = О (fxs + р); а также для уравнений с возмущающей силой от гауссова кольца, моделирующего влияние от второго тела: С (а) 1, CQ(O) ф 0 [p.s + fJ-p + fyx). Во всех случаях ро = 0, а начальные данные для третьего тела соответствовали оскулирующей круговой орбите с единичным радиусом.
Из рисунка видно, что учет вековых эффектов позволяет существенно повысить точность представления движения. Низкая точность вблизи а = 0.5, по-видимому, объясняется тем, что под сильным гравитационным влиянием от второго тела орбита тре Если третье тело — Юпитер, то это соответствует интервалу 1000 лет тьего значительно отличается от круговой, поэтому полученные нами оценки вековых эффектов становятся слишком грубыми. Кроме того, следует заметить, что несистематические всплески для а 0.3 главным образом обусловлены влиянием ошибок округления, которые, в свою очередь, вызваны значительным уменьшением шага интегрирования и потому увеличением объема вычисления, что требуется для интегрирования короткопериодических возмущений.
Расхождения решений для других уравнений значительно больше. Хотя для уравнений, где С (а) = 1 и С (а) 1 (гауссово кольцо), с уменьшением а точность их решений повышается, причем для С (а) 1 ошибка [Дх убывает быстрее. Тогда как интегрирование невозмущенных уравнений дает почти одно и то же довольно большое расхождение решений независимо от а.
Интересно отметить, что при а = 0.3 (Марс) моделирование влияния от планеты с помощью гауссова кольца дает тот же результат, что и простое включение массы планеты в центральное тело. Поэтому в данном случае не имеет никакого смысла использовать сложные уравнения с гауссовым кольцом вместо простых уравнений задачи двух тел, но с модифицированным гравитационным параметром.
К сожалению, оценки вековых эффектов (2.23), полученные нами для кругового случая, оказываются совершенно бесполезными даже для умеренно эксцентричных орбит. В то же время аналитические выкладки для 8 в общем случае весьма затруднительны. Далее мы предлагаем численный алгоритм для приближенного вычисления 5.
Идея алгоритма состоит в оценке 5 по отклонению решений, полученных из усредненных и исходных уравнений, на малом интервале времени порядка одного оборота. Затем на оставшемся, более длительном интервале интегрируются модифицированные усредненные уравнения, учитывающие вековые эффекты в угловой переменной, за счет чего и удается повысить быстродействие численного процесса.
Очевидно, проблема низкого быстродействия вследствие короткопериодических возмущений становится наиболее ощутимой главным образом в эволюционных задачах, когда исследование движения проводится на длительных временных интервалах порядка 100 оборотов объекта и более. Поэтому ценность предлагаемого алгоритма будет тем выше, чем длиннее исследуемый интервал времени. Впрочем, следует заметить, что хотя алгоритм и приближенный, и не позволяет полностью учесть вековые эффекты в быстрой переменной, тем не менее с его помощью все же удается повысить точность представления движения на несколько порядков.
Чтобы исследовать точность представления движения уравнениями (2.29), мы провели численный эксперимент, подобный описанному выше, но на примере пространственной задачи. Напомним, что проблема понижения быстродействия численного интегрирования особенно важна при моделировании спутниковой динамики, которая обременена короткопериодическими возмущениями от близких массивных спутников типа галилее-вых. Чтобы приблизить эксперимент к условиям спутниковой задачи, мы взяли начальные условия для третьего тела, соответствующие орбите одного из далеких спутников Юпитера, Гималии, с орбитальными элементами е = 0.16, г = 28 и приняли fis = 1) /хр = 10 4. Интегрирование выполнялось интегратором Эверхарта 15-го порядка с переменным шагом на интервале времени оборотов объекта (100 лет).
Важность упомянутой выше проблемы наглядно демонстрирует рисунок 2.4. Так, если мы интегрируем орбиту Гималии, только учет влияния от Ио будет требовать уменьшение шага интегрирования, а потому понижения быстродействия, в 8 раз.
На рисунке 2.5 приведены результаты эксперимента, отклонения решений упрощенных уравнений относительно истинного решения. Как показывает рисунок, качественно они такие же как и на рисунке 2.3 для планетной задачи. Хотя здесь для модифицированных уравнений величина 5 вычислялась по изложенному выше алгоритму, причем на интервале At, равном 1 обороту. Из рисунка видно, что учет таким способом остаточных вековых эффектов в быстрой переменной позволяет повысить точность представления движения на несколько порядков.
Мы также исследовали эффективность применения уравнений (2.29) для моделирования реального движения далекого спутника Гималии. В структуру возмущений мы включили влияние всех галилеевых спутников и Солнца. Причем координаты массивных спутников вычислялись по формулам кругового движения в плоскости экватора Юпитера13, а координаты Солнца — из эфемериды DE405 (Standish, 1998). В усредненных уравнениях влияние галилеевых спутников моделировалось с помощью гауссовых колец. На рисунке 2.6 показано поведение ошибок моделирования при использовании различных упрощенных уравнений. В частности, из рисунка видно, что модификация усредненных уравнений позволяет повысить их точность на несколько порядков.
Ограниченная задача трех тел. Астероидная задача
Чтобы детально исследовать возможности рассмотренных методов, мы применили их к решению упрощенной обратной задачи, где движение спутника моделировалось на основе формул задачи двух тел (приложение 2). Модель представляла движение спутника Адрастеи в пространстве угловых координат. В качестве орбитальных параметров были взяты компоненты начального вектора динамического состояния q = (х0,х0), предварительно полученные из наблюдений методом Лапласа (Escobal, 1965). При этом параметры q соответствовали почти круговой орбите с большой полуосью а = 8.68-Ю-4 а.е. и эксцентриситетом ё = 0.0161. На основе вектора q мы моделировали (точные) наблюдаемые положения р = (а,6) в моменты реальных наблюдений (N = 90), рассредоточенных на концах интервала времени 12 лет как две почти равновесные группы. Начальная эпоха to была взята в середине временного интервала.
Варьируя элементы а и е и получая различные начальные приближения qo, мы определяли орбитальные параметры q по р в соответствии со схемой Гаусса-Ньютона применительно к (4.22). Таким образом, выполнив многочисленный ряд экспериментов, мы обнаружили, что итерационный процесс сходится далеко не при всех q0, а если и сходится, то не всегда к истинному решению. Причем область сходимости к ожидаемым элементам due оказывается довольно малой. По-видимому, причина расходимости итерационного процесса (даже при хороших начальных приближениях q0, обеспечивающих малые значения целевой функции) главным образом связана здесь с сильной овражностью и сложной структурой гиперповерхности, задаваемой (q).
Для улучшения сходимости итерационного процесса мы применили демпфированный метод Гаусса-Ньютона (4.27). На рисунке 4.7 показаны результаты, полученные при h = 0.01. Здесь а — относительное отклонение большой полуоси а от истинного значения а, т.е. а — а(1 + а). Точками обозначены решения, к которым сходился итерационный процесс при различных начальных вариациях а и е (внутри прямоугольника с пунктирной границей). Числа слева от точек означают, сколько раз итерационный процесс сходился к соответствующему решению, тогда как числа справа представляют среднеквадратическую ошибку а = л/S/N, получаемую в этом решении. При этом рассматривалось 2000000 начальных приближений, из которых сошлось только 0.4%, а из них, в свою очередь, большая часть к решениям внутри области варьирования элементов.
Итак, в рассматриваемой задаче также имеет место множество решений, среди которых, однако, доставляет абсолютный минимум только соответствующее исходным орбитальным элементам а и ё. Следует заметить, что почти все минимумы вдоль а распределены равномерно с шагом Аа « 9.23 Ю-5. Нетрудно проверить, что это рас 106 пределение хорошо согласуется с оценкой (4.17) при четных т. Причина отсутствия минимумов для нечетных т состоит в особом выборе начальной эпохи. Вернемся к рисунку 4.2. Его нижний график показывает, что в круговой задаче при начальной эпохе, находящейся в середине временного интервала (как в нашем случае), решения, доставляющие минимум целевой функции, должны лежать на орбите в почти диаметрально противоположных точках. Варьируя только большую полуось и эксцентриситет, мы фактически исключили те начальные приближения, которые могли бы дать минимумы в противоположной части орбиты, соответствующие нечетным т в оценке (4.17). В то же время верхний график на рисунке 4.2 показывает, что выбор эпохи внутри одной из групп локализует минимумы, что естественно удобно для их поиска.
Далее, поскольку мы предполагаем, что моделируемые наблюдения безошибочны, в абсолютном минимуме о = 0. В соседних минимумах а = 0.23". Поэтому при наличии в наблюдениях случайных ошибок с дисперсией порядка этой среднеквадратической величины, минимум а для истинной орбиты может быть сравним с соседними аналогами, причем настолько, что среди них он не будет узнаваем как соответствующий истинной орбите.
Многократно выполнив процедуру определения орбиты при различных выборках нормально распределенных ошибок, вносимых в моделируемые наблюдения, мы сравнили значения минимумов для истинной орбиты и одного из соседних аналогов. На рисунке 4.8 показаны вероятностные плотности Р их разностей Д т в случаях трех дисперсий ошибок s = 0.23", 0.46", 0.67". Положительное значение разности означает, что минимум, соответствующий истинной орбите, меньше соседнего.
В частности, из рисунка видно, что при достаточно малых ошибках (s = 0.23") разности Асг принимают положительные значения, среди которых наиболее вероятны Д т = 0.02-0.03". То есть, таким образом, параметры, доставляющие наименьший минимум, могут рассматриваться как наилучшие оценки истинных орбитальных параметров. В отличии от этого случая в двух других имеют место отрицательные разности, вероятности появления которых сравнимы с вероятностями положительных разностей. Это как раз говорит о том, что по значениям двух близких минимумов невозможно с уверенностью оценивать степень правдоподобности соответствующих им параметров.
Установить статусы соседних минимумов удается, если имеется большое количество (причем достаточно точных) наблюдений. В противном случае это предприятие оказывается ненадежным. На том же рисунке 4.8 приведены результаты эксперимента для 20 моментов наблюдений (по 10 на каждую группу). Отсюда мы видим, что большая часть разностей во всех случаях (N = 20) концентрируется около нуля. При этом вероятность появления отрицательных разностей довольно высока.
Как отмечалось выше, сходимость итерационной схемы Гаусса-Ньютона может быть достигнута, если в нее ввести уменьшающий множитель h (4.27). Для получения численных результатов в предыдущем разделе мы использовали h = 0.01. Однако, учитывая низкий процент сходимости решений (0.4%), можно полагать, что для более уверенного
Демпфированный метод Гаусса-Ньютона
Естественно, коэффициент к можно также рассматривать как показатель степени нелинейности задачи. Очевидно, в линейном случае он равен нулю. Подобный показатель нелинейности был введен в работе (Черницов и др., 2006).
На рисунке 5.14 для каждого спутника приведено максимальное значение коэффициента хтах из вычисленных в 12 вершинах эллипсоида, построенного по соотвествую-щей ковариационной матрице при ка = 4.5. Как видно из рисунка, почти все именованные спутники (кроме Кале) имеют достаточно малые коэффициенты нелинейности, которые составляют величины порядка 0.001 и меньше. Размеры вероятностных областей
Максимальное значение коэффициента х (белая заливка). Серой заливкой показана характеристика нелинейности при исключении из рассмотрения двух противоположных вершин доверительного эллипсоида линеаризированной задачи с максимальными значениями х. для именованных спутников, как было показано, сравнительно небольшие и поэтому вариации областей за счет такой нелинейности можно считать незначительными.
В то же время значения показателя нелинейности хтах для некоторых неименованных спутников (открытых в 2003 г.) достаточно большие и поэтому для определения точности орбитальных параметров этих спутников требуется прибегать к нелинейному оцениванию. С этой целью были пересмотрены точностные характеристики для спутников J02, J03, J04, J09, J10, J12 и J14 группы S/2003, у которых хтах 0.01.
Интересно заметить, что для спутников с большим хтах нелинейность существенно проявляется только в двух противоположных (наиболее удаленных друг от друга) вершинах доверительного эллипсоида линеаризированной задачи. Как раз они задают направление, вдоль которого вытягивается доверительная область. На рисунке 5.14 серой заливкой представлены значения хтах при исключении из рассмотрения значений х в указанных (аномальных) вершинах эллипсоида, откуда видно, что в оставшихся вершинах влияние нелинейности незначительно, и показатель хтах не превышает 0.002.
Итак, в результате применения схемы (5.20) для построения вероятностных областей были получены нелинейные оценки точности орбитальных параметров для рассматриваемых спутников группы S/2003. На рисунке 5.15 в качестве примера приведены вероятностные области линейного и нелинейного оценивания для спутника S/2003 J04, ориентированные относительно собственных векторов w соответствующей ковариационной матрицы. Здесь даются только проекции вдоль собственного вектора W6
Несмотря на то, что разброс нелинейных оценок существенно больше вдоль соствен-ных векторов W1-W5, неопределенности в позиционных орбитальных элементах, полученные по этим оценкам, оказываются заметно меньше (рисунок 5.16). Тем не менее интервалы неопределенности в элементах остаются еще довольно большие и они имеют взаимные пересечения. Учитывая это, попутно мы попытались среди всех далеких спутников найти такие, возможные элементы которых на 1 января 2009 г. попадают в область неопределенности другого спутника, что несомненно представляет интерес с точки зрения идентификации объектов при их наблюдении в будущем23. Примечательно, что таких объектов не оказалось. Хотя заметим, что если бы для спутников использовались только линейные оценки, то обнаружилось бы пересечение областей неопределенности двух спутников группы Ананке, Ортозие и S/2003 J04.
Наконец, следует обратить внимание на то, что для спутников S/2003 J09 и S/2003 J14 линейные и нелинейные оценки (рисунок 5.16) хорошо согласуются, что не удивительно, так как именно эти спутники среди рассматриваемых имеют наименьшие коэффициенты нелинейности (рисунок 5.14).
Были также пересчитаны величины smax (5.24), характеризующие разброс орбит в пространстве измеряемых величин. Как видно из рисунка 5.17, для некоторых спутников они значительно меньше в сравнении с их линейными аналогами (рисунок 5.12), в особенности для S/2003 J04. Несмотря на это, ввиду все еще больших smax проблема потери объектов J02, J03, J04, J10 и J12 группы S/2003, так или иначе, остается.
При исследовании динамики новых спутников на длительных интервалах времени было обнаружено, что некоторые возможные орбиты спутника S/2003 J02 выходят за пределы гравитационной сферы Юпитера, иначе говоря, имеется вероятность, что объект станет астероидом. По предварительным линейным оценкам вероятность того, что уже за 100 лет спутник сменит свое амплуа, составляет приблизительно 0.06 (рисунок 5.18) (Авдюшев, 2008). Однако уже по нелинейным оценкам эта вероятность оказывается на порядок меньше и составляет 0.005.
Впрочем, следует иметь в виду, что представленные выше результаты основаны на оценках орбитальных параметров, хорошо согласующихся с полученными другими авторами (Sheppard et al., 2003; Emelyanov, 2005) (см. также приложение 4), в то время как существуют альтернативные оценки (XQ HXJB приложении 6). Их появление фактически связано с известной проблемой двузначности в предварительном определении орбит (Субботин, 1968), которая формально вызвана наличием двух решений уравнения
Линейные и нелинейные оценки неопределенностей в позиционных элементах спутников группы S/2003 (наклонение отнесено к геоэкватору)
Лагранжа, в частности, при определении спутниковых орбит (Charlier, 1910; Crawford, 1913). Причем обе получаемые орбиты одинаково хорошо представляют наблюдения.
При решении обратной задачи с использованием большой выборки наблюдений разумным критерием выбора нужного решения, соответствующего реальной орбите, является малость среднеквадратической ошибки. Однако, если в обоих случаях средне-квадратические ошибки сравнимы, выбор решения без привлечения дополнительной информации об орбите оказывается под вопросом. Такая ситуация вполне возможна для объектов, наблюдавшихся на короткой дуге, поэтому для поиска альтернативных решений мы ограничились спутниками семейства S/2003.
Нам удалось найти лишь 12 альтернативных орбит (таблица 5.1), из которых 7 оказались гиперболическими, т.е. согласно новым оценкам уже 7 объектов могут являться астероидами, временно захваченными Юпитером.
