Содержание к диссертации
Введение
1 Астероидно—кометная опасность 18
1.1 Осознание проблемы 18
1.1.1 Падение на Землю космических тел 19
1.1.2 Возможность глобальной катастрофы 20
1.1.3 Анализ риска 21
1.2 Объекты, сближающиеся с Землей 23
1.2.1 Минимальное межорбитальное расстояние 23
1.2.2 Классификация объектов, сближающихся с Землей 26
1.3 Туринская и Палермская шкалы 31
1.3.1 Туринская шкала 31
1.3.2 Палермская шкала 32
2 Определение орбиты малого тела Солнечной системы 35
2.1 Метод перебора плоскостей орбиты малого тела 38
2.2 Определение орбиты по двум гелиоцентрическим положениям и моментам времени методом Гаусса 40
2.3 Дифференциальный метод улучшения орбиты 43
2.4 Вычисление матрицы коэффициентов условных уравнений 49
2.5 Результаты 51
3 Методы оценки вероятности столкновения малого тела с Землей 54
3.1 Метод Монте-Карло 56
3.1.1 Основная идея 56
3.1.2 Реализация нормального закона распределения 56
3.1.3 Погрешность метода 58
3.1.4 Преимущества и недостатки 59
3.2 Линейные методы оценки вероятности 60
3.2.1 Описание метода 60
3.2.2 Учет корреляции с компонентами скорости 61
3.2.3 Учет гравитационного притяжения Земли 62
3.2.4 Спектральное разложение матрицы 64
3.2.5 Преимущества и недостатки 65
3.3 Метод плоскости цели 66
3.3.1 Плоскость цели 66
3.3.2 Описание метода 68
3.3.3 Преимущества и недостатки 69
3.4 Метод вариации одного параметра и метод LOV 69
3.4.1 Основная идея 69
3.4.2 Описание метода 70
3.4.3 Преимущества и недостатки 73
4 Новый линейный метод оценки вероятности 76
4.1 Система координат 78
4.1.1 Случай круговой орбиты 78
4.1.2 Случай эллиптической орбиты 79
4.1.3 Переход из (,г),М) в эклиптическую систему координат 80
4.1.4 Переход из эклиптической системы координат в криволинейную систему 82
4.2 Вычисление вероятности 87
4.2.1 Определение нормальной матрицы 87
4.2.2 Учет корреляции со скоростями 89
4.2.3 Нахождение области О в криволинейной системе координат 90
4.2.4 Спектральное разложение матрицы 91
4.2.5 Вычисление интегралов 92
4.3 Результаты 94
4.3.1 Вероятности, вычисленные методом Монте-Карло 99
4.3.2 Сравнение с линейными методами оценки вероятности 101
4.3.3 Сравнение с методом вариации среднего движения и LOV 104
4.3.4 Случай астероида 2010 RF12 107
Заключение 110
Список литературы 113
- Возможность глобальной катастрофы
- Определение орбиты по двум гелиоцентрическим положениям и моментам времени методом Гаусса
- Реализация нормального закона распределения
- Переход из эклиптической системы координат в криволинейную систему
Введение к работе
Актуальность темы. На данный момент интерес к проблеме астероидно-кометной опасности достаточно высок, особенно это связано с падением на Землю Чебаркульского метеорита. Данное явление показало, что даже не очень большие по космическим масштабам астероиды могут представлять определенную опасность для инфраструктуры и жизни людей, в случае столкновения таких объектов с Землей. Астероидно-кометная опасность выделяется из круга задач противодействия различным катастрофам, поскольку столкновение с большим космическим телом может привести к гибели всего живого на планете.
Для предотвращения этой угрозы необходимо иметь комплекс методов для определения орбит малых тел Солнечной системы и степени их опасности для Земли. Оценка степени опасности того или иного объекта напрямую связана с задачей определения вероятности столкновения объекта с Землей. В связи с возросшим темпом получения наблюдений астероидов, сближающихся с Землей, необходимо иметь быструю (линейную) и надежную методику для оценки вероятности столкновения малых тел с Землей, поскольку после каждого нового наблюдения объекта значение вероятности меняется и ее приходится оценивать заново. Существующие линейные методы оценивают вероятность столкновения с недостаточной точностью в случае, когда потенциальное столкновение происходит далеко от номинального положения астероида.
Цель и задачи работы. Целью данной работы являлась разработка быстрой методики обнаружения потенциальных столкновений с Землей вновь открытых малых тел Солнечной системы и вычисление вероятности этих событий, которую можно использовать для случая, когда моменты потенциальных столкновений далеки от эпохи наблюдения. Для этого необходимо решить следующие задачи:
Разработать линейный метод оценки вероятности столкновения астероида с Землей, который можно использовать для случая, когда потенциальное столкновение далеко от номинального положения астероида.
Произвести сравнение всех используемых методов оценки вероятно
сти столкновения астероидов с Землей.
Научная новизна работы.
Разработана методика определения орбит вновь открытых малых тел Солнечной системы, основанная на переборе плоскостей орбитального движения.
Разработана оригинальная криволинейная система координат, связанная с орбитой малого тела.
Решена задача нахождения уравнения формы Земли в новой разработанной системе координат.
Разработан новый линейный метод оценки вероятности столкновения малого тела с Землей, использующий оригинальную криволинейную систему координат.
Проведено детальное сравнение используемых в мире методов оценки вероятности столкновения. Показано, что в определенных случаях нелинейные методы вариации одного параметра дают менее точные оценки вероятности столкновения, чем линейные методы.
Практическая ценность. Важной задачей является оперативное определение орбиты и вероятности столкновения малого тела Солнечной системы с Землей. Актуальной научной и практической задачей является создание и поддержка списка угрожающих Земле небесных тел.
Создан комплекс программ, позволяющий определять орбиты вновь открываемых малых тел Солнечной системы и оперативно оценивать вероятности их столкновения с Землей. Данный комплекс применяется в работах, связанных с астероидно-кометной опасностью, в лаборатории малых тел Солнечной системы ИПА РАН.
Публикации по теме диссертации. Основные результаты по теме диссертации опубликованы в 16 работах. Во всех работах диссертант участвовал в постановке задачи, им были решены поставленные задачи и проведен
анализ полученных результатов. Перечень работ приведен в конце автореферата.
Результаты, выносимые на защиту.
Разработан новый линейный метод оценки вероятности столкновения малого тела Солнечной системы с Землей, использующий криволинейную систему координат, связанную с номинальной орбитой малого тела.
Результаты сравнения методов показали преимущество разработанного линейного метода оценки вероятности столкновения по сравнению с линейным методом, использующим декартову систему координат, и с методом плоскости цели.
Теоретически и практически показано, что в определенных случаях нелинейные методы вариации одного параметра дают менее точные оценки вероятности столкновения, чем линейные методы.
Статистически получено, что метод Line Of Variations Sampling является более надежным, чем метод вариации среднего движения.
Апробация работы. Основные результаты диссертации докладывались на семинарах ИПА РАН, ГАО РАН, ИНАСАН и на семинарах кафедры небесной механики СПбГУ. Результаты также докладывались на конференциях: "IV Пулковской молодежной астрономической конференции" (ГАО РАН, 2012), Всероссийской астрометрической конференции "Пулково-2012" (ГАО РАН, 2012), молодежной конференции "Фундаментальные и прикладные космические исследования" (ИКИ РАН, 2013), и на международных конференциях "Околоземая астрономия-2013" (Туапсе, 2013), "Asteroids, Comets, Meteors" (Хельсинки, 2014), "Journees" (Санкт-Петербург, 2014), "GAIA-FUN-SSO 3" (Париж, 2014), "Planetary Defence" (Рим, 2015).
Структура диссертации. Диссертация состоит из четырех глав, введения, заключения и списка литературы. Она изложена на 118 страницах
Возможность глобальной катастрофы
Тунгусский феномен наглядно показывает, что встреча Земли с телом размером в сотню или более метров может обернуться региональной катастрофой. Однако о том, что столкновение Земли с достаточно крупным космическим объектом грозит человечеству глобальной катастрофой, стало известно только после того, как в разных странах было осуществлено математическое моделирование последствий тотального ядерного конфликта [20]. Как оказалось, наиболее вероятным последствием такого конфликта будет эффект долговременной "ядерной зимы", возникающий вследствие выброса в атмосферу колоссального количества мелкодисперсной пыли и пепла. Пыль и пепел, достаточно быстро распространятся во всей атмосфере Земли или ее значительной части, где могут оставаться в течение многих месяцев и даже лет, экранируя солнечный свет. Произойдет глобальное нарушение динамики всей атмосферы Земли. Температура верхних слоев атмосферы резко возрастет, при это температура у поверхности Земли упадет на десятки градусов. Растения прекратят синтез органического вещества. Теплокровные животные погибнут от снижения температуры и недостатка пищи. Почва, внутренние водоемы и поверхностные слои морей и океанов окажутся отравленными кислотными дождями. Спустя месяцы, после очищения атмосферы и восстановления ее циркуляции возникнет парниковый эффект вследствие увеличения содержания в атмосфере углекислого газа. В результате может произойти глобальное самоподдерживающееся повышение температуры Земли на несколько градусов. Температурные колебания могут иметь драматические последствия для всей окружающей среды и обернуться гибелью большинства обитателей Земли.
Причина описанной глобальной катастрофы практически никак не связана с ядерной природой взрыва. Нечто подобное должно иметь место при ударе о Землю тела размером в 1-2 км и более. Из образовавшегося кратера в десятки километров в диаметре выбрасывается большое количество вещества, в том числе в виде частиц микронного размера. Огненные бури, колоссальные лесные пожары способны охватить целые регионы. Нарушения в динамике атмосферы и инсоляции могут привести к массовой гибели обитателей Земли.
В работе Чапмана и Моррисона 1994 года [30] была оценена средняя вероятность смерти человека от удара космического тела в течении жизни. В таблице 1 приведены вероятности смерти людей в США по различным причинам. Отметим, что очень сложно оценить вероятность смерти по причине падения метеорита на Землю с неплохой точностью, поэтому в таблице приведены верхний и нижний пределы вероятности, а также средняя оценка этого значения.
Стоит также отметить, что эти оценки делались на основе выражения зависимости числа жертв N от кинетической энергии ударника Е , выраженной в мегатоннах: N = 2 1073. Коэффициент перед Ek зависит от средней плотности населения Земли на км . В 1994 году численность населения Земли была равна 5.6 миллиарда человек и средняя плотность населения 10 человек/км . В настоящий Таблица 1: Вероятность смерти по различным причинам (в США)
Вычисленные указанным способом шансы умереть от удара космического тела в среднем такие же, как от аварии самолета или от наводнения, и значительно ниже, чем вероятность погибнуть при пожаре или в автомобильной катастрофе. Более того, пока неизвестны случаи, чтобы кто-то погиб от удара космического тела (из рассказов очевидцев Тунгусского события можно заключить, что были смерти лишь из-за сильного психологического действия взрыва). Но можно заметить, что большое количество погибших в одном событии, например при аварии авиалайнера, обычно воспринимается как большее несчастье, чем постепенная и частая гибель людей, например в автоавариях. Гибель около 3000 людей из-за нападения террористов 11 сентября 2001 г. имела гораздо более серьезные последствия (включая экономические, политические и международные), чем то же число погибших в авариях на шоссе. В этом отношении удары космических тел занимают особое положение, так как могут быстро убить огромное количество населения. Различные наиболее разрушительные природные катастрофы, случившиеся в истории, несмотря на их мощь, убили лишь относительно небольшую часть всего населения: землетрясения — 2 млн человек, циклоны — 300000, обвалы и оползни — 100000, цунами — 100000, вулканы — 90000, лавины — 20000. Удар космического тела может убить всех.
При анализе объектов, угрожающих столкновением с Землей, одной из наиболее важных характеристик является MOID (в англоязычной литературе эта аббревиатура расшифровывается как Minimum Orbital Intersection Distance). Эта величина представляет собой минимальное расстояние между двумя кеплеровыми орбитами небесных тел, не принимая во внимание положения объектов на этих орбитах. Как таковой MOID может рассматриваться как первоначальный индикатор опасности столкновения небесного тела с Землей. Большой MOID между орбитами астероида и Земли означает, что астероид не столкнется с Землей в ближайшем будущем, в то время как за астероидами, орбиты которых имеют маленький MOID, нужно тщательно следить, поскольку они могут иметь тесные сближения или столкновения с Землей. Впервые этот термин был введен в работе Боуэлла и Мюйнонена [29].
Давайте рассмотрим две кеплеровы орбиты. Как хорошо известно эти орбиты являются коническими сечениями, которые могут быть ограниченными (эллипс и его частные случаи круг и отрезок), или неограниченными (парабола, гипербола и прямая). Более того, эти два конических сечения имеют общий фокус в центре Солнца.
Пусть Е\ и Е два множества из 5-ти элементов, которые описывают выбранные кеплеровы орбиты. В качестве этих пяти параметров могут ис пользоваться различные системы элементов орбит. Одним из возможных вариантов может выступать Ек = {а , е&, І&, Г , }, (к = 1,2), где ( — большая полуось орбиты, е& — эксцентриситет, і к — угол наклона, Qk — долгота восходящего узла и Шк — аргумент перицентра. Обозначим за V\ и г 2 элемент, отвечающий за положение объекта на кеплеровой орбите. Это может быть, например, истинная, эксцентрическая или средняя аномалии. Если обозначить за ф\ = 0і(і,г і), ф = 02( 2,) декартовы координаты первого и второго тела на орбитах с компонентами (жі, yi, Z\) и (х2-, у 2-, 2), соответственно, то величина MOID в общем случае может быть определена следующим образом:
Так как модуль — функция всегда положительная, тогда инфимум, являясь максимальной нижней границей множества, существует всегда, следовательно и MOID по определению существует всегда. Если инфимум достигается на самих множествах, т.е. в нашем случае существуют по одной точке на выбранных орбитах, расстояние между которыми равно MOID-y их орбит, в таком случае MOID является абсолютным минимумом Евкли-дового расстояния между точкой на первой орбите и точкой на второй. Для случая двух ограниченных орбит MOID всегда будет достигаться, поскольку множества точек пространства, которые они представляют, компактны, и следовательно, будет существовать минимум функции \/{х\ — Х2)1 + (у\ — У2)1 + {z\ — Z2)2, хотя он может быть и не единственным. Например, два эллипса, лежащие в разных плоскостях, могут иметь две точки, в которых достигается минимальное расстояние. В частности, они могут пересекаться в двух точках, что будет означать, что MOID этих орбит равен нулю. Для случая двух окружностей, лежащих в одной плоскости и имеющих общий центр, имеется континуум пар точек, в которых достигается минимальное расстояние между орбитами.
Определение орбиты по двум гелиоцентрическим положениям и моментам времени методом Гаусса
Как можно заметить, задача нахождения орбиты малого тела так или иначе связана с задачей определения топоцентрических расстояний объекта. В методах Гаусса и Лапласа, а также в их модификациях, нахождение топоцентрических расстояний производится с использованием метода последовательных приближений. При этом могут возникнуть проблемы со сходимостью итераций. Основная цель этого метода — избежать применения последовательных приближений для нахождения топоцентрических расстояний. Вместо нахождения невозмущенной кеплеровой орбиты малого тела, проходящей через три наблюдения, мы находим первое приближение орбиты, а затем, улучшаем орбиту дифференциальным методом.
Заметим, что топоцентрические расстояния - непосредственно зависят только от угла наклона орбиты і, долготы восходящего узла Q и от наблюдений lj, которые нам известны. Для наглядности запишем уравнение (2.1) в следующем виде:
В предложенном методе начальное приближение орбиты малого тела находится перебором плоскостей орбитального движения, т.е. перебором угла наклона і и долготы восходящего узла Q плоскости орбиты [6]. Прежде всего, выбираются 2 опорных наблюдения (обычно это первое и последнее наблюдение, но возможны и другие комбинации, в случае, если известно, что ошибки некоторых наблюдений отличаются от остальных). Затем для каждой рассматриваемой плоскости выполняется следующее. Значения то-поцентрических расстояний pj находятся как длины отрезков векторов направления lj на объект от наблюдателя до его пересечения с плоскостью, с выбранным наклоном и долготой восходящего узла. Величины pj определяются формулой: где N = (sin і sin Г2, —sin і cosil, cos і) — вектор нормали к плоскости орбиты.
Следует отметить, что близкий по идеологии подход применялся в работе Виртанена и Мюйнена [52]. Ими в качестве перебираемых параметров были выбраны геоцентрические расстояния для двух моментов времени.
Далее, на моменты времени наблюдений вводятся поправки за аберрацию (t A = tj — -pj) , где с — скорость света. Затем, используя метод Гаусса для определения элементов орбиты по двум гелиоцентрическим положениям и моментам времени [22], находится орбита. Наконец, вычисляются разности между наблюденным и вычисленным положением тела (О — С), и определяется значение среднеквадратической ошибки представления наблюдений: п
Величина о" показывает насколько хорошо эта орбита описывает данный набор наблюдений. Если каким-либо образом известны веса наблюдений или они предполагаются известными, то их можно учесть уже непосредственно в выражении для а [11]. Стоит также отметить, что тут была выбрана Евклидова норма для нахождения т, хотя возможно использование любой другой нормы векторного конечномерного пространства, например Гёльдерова норма ж = ( ж р)р, частным случаем которой яв i
Среди всех вычисленных значений среднеквадратической ошибки а вы бираем наименьшее. Элементы орбиты, дающие наименьшее т, берутся за начальное приближение орбиты малого тела. В случаях, когда имеется несколько орбит с близкими значениями о" и сильно различающимися наборами элементов, рассматриваются несколько возможных вариантов орбиты.
В том случае, если полученные элементы орбиты не улучшаются дифференциальным методом улучшения орбит, производится поиск более точной невозмущенной орбиты. Для этого шаги по значениям углов наклона и долготы восходящего узла уменьшаются, и рассматриваются другие варианты выбора опорных наблюдений, отличных от крайних.
Главное преимущество описанного подхода — всегда наличие некоторого начального приближения. Классические методы, в этом плане, являются не очень надежными. В случае возникновения проблем, например, со сходимостью итераций, на выходе мы не получаем ничего. Второе преимущество данного метода, что если орбита не может быть улучшена, мы можем рассматривать множество орбит, которые удовлетворяют наблюдениям с заданной точностью.
Пусть имеется 2 гелиоцентрических положения тела Г\ = (xi,yi,Zi) и г2 = ( 2,2/25) в моменты времени t\ и 2, соответственно. Задача нахождения орбиты в первую очередь заключается в том, чтобы найти параметр орбиты р [22]. Как будет показано далее, этого достаточно, чтобы определить оставшиеся элементы орбиты. Впервые задачу нахождения параметра орбиты по двум гелиоцентрическим положениям и моментам времени решил Гаусс в 1809 г.
Реализация нормального закона распределения
Как уже отмечалось выше, оставшийся трехмерный интеграл вычисляется численно. Однако для того, чтобы вычислить его с неплохой точностью потребуется много машинного времени. Для существенного уменьшения времени вычисления существует следующий подход. Используется спектральное разложение матрицы N = U N U, где N — диагональная матрица, на диагонали которой стоят собственные числа матрицы N, и U — ортогональная матрица, столбцы которой являются ортонормированны-ми собственными векторами матрицы N. Как это делалось ранее, переходим в новое пространство z = Ut/, где у — вектор декартовых координат. Поскольку матрица U является ортогональной, то эту матрицу можно рассматривать как матрицу поворота. Таким образом, это преобразование является просто переходом посредством поворота в другую декартову систему координат, в которой координаты являются некоррелированными. В этом смысле любая декартова система координат одинаковым образом описывает функцию распределения возможных положений астероида.
После того, как мы перешли в новую систему координат, функция распределения ошибок положений астероида разбивается на произведение трех функций вида:
Поскольку произошел только поворот системы координат, то область, которую занимает Земля, осталась шаром с тем же радиусом Ле. Так как функция, от которой берется интеграл, всегда положительна, то если заменить область Земли на куб, с тем же центром и стороной 2Ле, то значение интеграла может только увеличиться. Выполним такую замену, взяв куб со сторонами, параллельными осям новой системы. Оставшийся трехмерный интеграл запишется в виде произведения: о zi+R(S 2 где (zf, zi, z\) — координаты центра Земли в новой системе координат, — элемент матрицы N в і-ой строке и і-ом столбце, Y\ — обозначение произведения. Вычисление численно трех однократных интегралов значительно отличается по времени от вычисления тройного. И чем больше требуется точность вычисления, тем больше будет различие во времени.
По причине замены области, по которой берется интеграл, на большую, значение вероятности, таким образом, было увеличено. Однако в связи с тем, что на практике ошибки (значения az.) двух из координат на один или несколько порядков меньше радиуса Земли, то если рассматривается момент, когда вероятность достигает своего максимального значения, это увеличение не является существенным для оценки вероятности столкновения.
Остается вопрос каким образом нужно просуммировать вероятность, чтобы получить полную вероятность столкновения, а не на конкретные даты. Как оказалось на практике, максимальное значение вероятности в окрестности некоторой даты будет очень близко к полной вероятности в окрестности этой даты и нет необходимости производить какое либо суммирование.
Основное преимущество линейных методов это скорость вычисления. При вычислении вероятности требуется проинтегрировать орбиту малого тела только один раз, в отличии от нелинейных методов, где требуется проинтегрировать больше нескольких тысяч орбит виртуальных астероидов.
Какие ограничения имеет данный метод? Во-первых, поскольку предполагается, что нормальный закон распределения ошибок сохраняется на всем рассматриваемом интервале времени, то это значит, что не должно быть тесных сближений с массивными телами, которые могут сильно отклонить закон распределения ошибок от нормального. Однако это не единственный источник нелинейности. Как можно заметить, предполагая нормальный закон в декартовой системе координат, получаемая область виртуальных астероидов будет эллипсоид в декартовой системе координат. Однако различия в среднем движении у виртуальных астероидов будут приводить к тому, что облако виртуальных астероидов будет со временем вытягиваться вдоль номинальной орбиты объекта. Если время, прошедшее после последнего наблюдения не очень большое и ошибки являются достаточно маленькими, то это облако не будет значительно отличаться от эллипсоида в декартовых координатах. Напротив, для случая больших ошибок или значительного времени от момента последнего наблюдения облако виртуальных астероидов сильно вытянется вдоль орбиты, и различия станут существенными. Этот факт делает ненадежным линейные методы, использующие декартову систему координат, для оценки вероятности столкновения небесных тел с Землей на больших интервалах времени.
Метод плоскости цели является одним из линейных методов оценки вероятности столкновения малого тела с Землей, однако он настолько распространен, что было принято решение выделить описание этого метода в отдельный пункт. Начнем с определения самой плоскости цели.
Плоскость цели или Ь-плоскость. Это классическая плоскость цели, которая использовалась в астродинамике с 60х годов прошлого века. Это плоскость, которая содержит центр Земли и перпендикулярна входящей асимптоте гиперболы, по которой движется тело внутри сферы действия Земли. Иными словами, нормаль b-плоскости сонаправлена с невозмущенной относительной скоростью астероида г оо до входа в сферу действия. Название Ь-плоскость пошло от так называемого b параметра (прицельного расстояния) — расстояния между асимптотой гиперболы и центром Земли, т.е. минимального расстояния от невозмущенной орбиты астероида и Землей.
Модифицированная плоскость цели (МПЦ). МПЦ это модифицированная плоскость в том смысле, что она перпендикулярна геоцентрической скорости астероида в момент самого тесного сближения вдоль номинальной траектории движения. Вместо того, чтобы отмечать пересечения невозмущенной орбиты, как это делалось в b-плоскости, мы отмечаем пересечения возмущенной траектории астероида с плоскостью цели, когда используем модифицированную плоскость цели. Термин МПЦ был введен впервые в 1999 году в работе Милани и Вальсеччи [44], хотя концепция использовалась раньше (например в работе Чодаса и Иоманса 1994 года [34]).
Главное отличие между двумя типами плоскости цели состоит в том, что отклонение астероида гравитационным полем планеты или гравитационное фокусирование напрямую связано с использованием модифицированной плоскости цели, в то время как при использовании b-плоскости этот факт скрывается. Для сближений на высоких скоростях или на больших расстояниях, отклонение слишком мало, так что различием между двумя плоскостями цели можно пренебречь.
При использовании МПЦ возникает трудность, так как близкие траектории отклоняются по разному. Когда отклонения значительны, это приводит к сильной нелинейности при проектировании положения астероида до сближения на плоскость цели. Более того, сама по себе ориентация МПЦ может значительно меняться при слабых изменениях траектории астероида. Этот факт приводит к еще одному источнику нелинейности. Конечно, в случае, когда происходит частичный захват астероида, асимптоты не существует и поэтому подход Ь-плоскости не может быть применен. В этом случае стоит использовать модифицированную плоскость цели.
Рассуждения выше позволяют нам оценить применимость Ь-плоскости и модифицированной плоскости цели. При очень слабом отклонении оба подхода, по существу, не различимы, в то время как для умеренных отклонений предпочтение той или иной плоскости стоит отдавать в зависимости от обстоятельств и требуемых целей.
Переход из эклиптической системы координат в криволинейную систему
Для вычисления вероятности столкновения астероида с Землей в момент времени t линейным методом, использующим введенную криволинейную систему координат, нужно определить нормальную матрицу N M ВО введенной системе в момент времени t. Пусть после определения орбиты малого тела у нас имеется вектор эклиптических координат и скоростей (жо,2/о, ZQ, ХО, уо, ZQ) на эпоху to и нормальная матрица N yz. Численным методом интегрируются уравнения движения в вариациях малого тела от эпохи to до момента времени t, и находится матрица изохронных произ водных Ф(о,) декартовых координат и скоростей. (В данной работе интегрирование проводилось методом Эверхарта 15-го порядка. Учитывались гравитационные возмущения от больших планет, Плутона, релятивистские эффекты от Солнца. Также возмущения от Земли и Луны рассматривались раздельно.) Матрица Ф(іо,і) записывается следующим образом: где t{ — время на і-ом шаге интегрирования (i = (1,&)). Такое разложение позволяет избежать возрастания ошибок в процессе интегрирования в следствии увеличения со временем коэффициентов вида -к?-. Однако для сохранения точности выполнять произведение матриц стоит, используя 32 знака в записи числа после запятой. Иначе эффект накопления ошибок может сильно исказить нормальную матрицу Nxyz вплоть до того, что ее определитель может стать отрицательным.
Для нахождения нормальной матрицы N M выполняется следующее. Сначала фиксируются 5 параметров оскулирующего эллипса. Вычисляются координаты астероида (,Г), М) в момент времени t во введенной системе координат как описано в предыдущем пункте ( = 7 = 0 поскольку сам астероид лежит на оскулирующем эллипсе своей орбиты). Скоростные компоненты , ), М находятся из соотношения: Г]
Поскольку значения скоростей ,?),М могут принимать любые значения от —оо до +оо, производится аналитически интегрирование по трем скоростным компонентам, как это описывалось в предыдущей главе в пункте про линейные методы оценки вероятности. При каждом интегрировании по скоростной координате элементы нормальной матрицы изменяются по следующим правилам: _ TljiUkl пи где rijk — элемент нормальной матрицы в j-ой строке и k-ом столбце, а / — номер координаты, по которой происходит интегрирование. Константа перед интегралом увеличивается в 4= раз.
Таким образом, шестимерная нормальная матрица N M становится трехмерной матрицей N1 м и шестимерная область G становится трехмерной. Предполагается, что в декартовой системе координат область G является шаром с радиусом Ле.
Однако, так же как и в линейном методе оценки вероятности в декартовой системе координат часто возникают ситуации, когда потенциальное столкновение происходит далеко от номинального положения астероида. В этом случае, интегрируя уравнения в вариациях, невозможно учесть гравитационную фокусировку виртуальных астероидов Землей, и радиус Земли стоит увеличить в Wl + y- раз, где vs — параболическая скорость у поверхности Земли (= 11.18 км/с), VQQ — скорость астероида при входе в сферу действия Земли.
Если номинальное положение астероида находится внутри сферы действия Земли в момент потенциального столкновения/:, возмущения от Земли уже не являются линейными и мы не сможем интегрируя уравнения движения астероида в вариациях найти функцию распределения. В этом случае стоит придерживаться следующей схемы. Интегрировать уравнения движения астероида в вариациях до входа в сферу действия Земли, а далее, интегрировать их без учета притяжения Земли или рассматривать невозмущенную гелиоцентрическую орбиту. Тогда радиус Земли стоит увеличить в л/1 + f- раз, для учета гравитационной фокусировки.
Пусть ( е,?7е,Ме) — координаты центра Земли в криволинейной системе. Тогда проекция области G на плоскость (,7?) есть круг с центром в (е, Г]е) и радиусом Ле. Заметим, что Лф С 1, а потому изменение М в этой области можно считать линейным. Это значит, что если ( декартова координата с центром в центре координат (,??), такая что е перпендикулярен и е и еч, то изменение координаты М можно выразить как изменение , умноженное на константу,
Константу d можно найти из следующих соображений. Если V модуль скорости астероида в эклиптической системе координат в точке, соответствующей Ме, то тг есть время, за которое астероид пройдет расстояние А( в этой области. Умножая это время на М (скорость изменения М) получаем изменение средней аномалии М. Поэтому константа будет равна а у .
Преобразуем систему координат ( ,ту,М) так, чтобы в новой системе область G имела бы форму шара с радиусом Ле. Новая система координат будет (,?7, -г). Чтобы определить нормальную матрицу Nf м в новой системе координат нужно умножить третий столбец и третью строчку матрицы N1 м на d. В этом случае элемент в третьем столбце и третье строке умножается на d . Константа перед интегралом тоже должна, таким образом, умножится на d.
Поскольку мы преобразовали нашу систему так, чтобы область Земли имела форму шара, можно воспользоваться подходом, примененным в линейном методе вычисления вероятности в декартовых координатах для увеличения скорости вычисления интеграла. Рассмотрим спектральное разложение матрицы Nt м = U N U. Перейдем к системе координат ( J7? d ) котРая является произведением матрицы U на (,77, -j). Т.к. матрица U является ортогональной, то область G остается шаром с тем же радиусом, при этом корреляции между элементами [с, ,Т) и—г J отсутствуют. Заменяя область Земли на куб с центром в центре Земли и стороной 2ЛФ, оставшийся тройной интеграл разбивается на произведение однократных вида: