A method for the fast MOID computation for two confocal heliocentric orbits

Full Text

Проблема астероидной опасности включает в себя множество различных задач, одной из которых является отнесение астероида к классу потенциально опасных [1, 2]. Согласно определению, используемому Лабораторией реактивного движения NASA (http://neo.jpl.nasa.gov/neo/groups.html), потенциально опасным астероидом (PHA - Potentially Hazardous Asteroid) считается такой, для которого одновременно выполняются два условия: орбита астероида имеет сближение с орбитой Земли на расстояние менее 0.05 а.е. (астрономических единиц) и диаметр астероида больше 150 метров (т.е. абсолютная звездная величина H 22). Для оценки расстояния между пересекающимися орбитами рассчитывается величина минимального расстояния между орбитами, сокращенно MOID (Minimum Orbital Intersection Distance). В настоящее время астероидов, принадлежащих к группам Аполлона, Амура и Атона, насчитывается более 12 000. В основном потенциально опасные астероиды принадлежат именно к этим группам [3-5]. По данным, представленным Лабораторией реактивного движения NASA (http://neo.jpl. nasa.gov/orbits) и Центром малых планет (http://www.minorplanetcenter. net/iau/Dangerous.html), на сегодняшний день насчитывается более 1 500 потенциально опасных астероидов. Таким образом, при получении новых данных об элементах орбит астероидов необходимо произвести расчеты для выявления потенциально опасных. В силу того, что расчеты необходимо проводить для значительного числа астероидов, время расчетов становится важным критерием при выборе метода оценки MOID. Кроме того, важна и точность получаемых оценок. Методы для оценки параметра MOID подразделяются на три группы: аналитические, численные и численно-аналитические. Среди аналитических методов классическими считаются методы, созданные К. В. Холшевниковым [6, 7] и G. F. Gronchi [8-10]. В работе К. В. Холшевникова [6] задача отыскания MOID сведена к решению тригонометрического уравнения восьмой степени, а также показано, что в общем случае дальнейшее упрощение задачи невозможно. В основе алгоритма, предложенного G. F. Gronchi, лежит быстрое преобразование Фурье. Задача сводится к алгебраическому полиному 16-й степени, действительные корни которого затем используются для отыскания критических точек функции, описывающей расстояние между орбитами [9, 10]. Как метод К. В. Холшевникова, так и метод G. F. Gronchi в результате реализации предоставляют информацию не только о точках минимума функции расстояния между двумя орбитами, но и о других критических точках этой функции. Методы, принадлежащие к группе численных, представляют собой итерационные алгоритмы. Задача нахождения MOID в них решается последовательным уточнением значения MOID. Достоинством численных методов является простота реализации и возможность настройки желаемой точно145 Д е р е в я н к а А. Е. сти и скорости работы. Схема работы методов этой группы состоит в следующем: сначала производится расчет расстояний между двумя точками на орбитах исследуемых тел, а потом проводится анализ дискретного представления функции расстояния между орбитами. Большинство современных методов являются гибридными, удачно сочетая в себе достоинства численных и аналитических методов [9, 11-15]. Численно-аналитические методы предоставляют высокую скорость с сохранением точности расчетов. Их общую структуру можно описать следующим образом. Аналитическими методами задача сводится к упрощенному виду, чтобы в дальнейшем применить численные методы для получения результата. В работах [1, 12, 15-20] используется именно такой подход. В данной работе представлен метод, позволяющий производить быструю оценку параметра MOID для двух тел на конфокальных эллиптических орбитах. Метод относится к группе численно-аналитических методов и предназначен для использования в случаях, когда требуется провести расчет параметра MOID для значительного числа астероидов, затратив при этом минимум времени. К примеру, такие расчеты необходимы в задаче оценки потенциальной опасности небесного объекта при оценке вероятности столкновения астероидов с Землей. Оценка параметра MOID на основе данных наблюдений астероидов проводится для первичного отбора потенциально опасных астероидов и выделения астероидов, чьи орбиты расположены близко к орбите Земли. Для отобранных астероидов проводится более тщательное интегрирование уравнений движения. Кроме того, выделение неопасных астероидов необходимо по той причине, что процедура оценки вероятности столкновения астероида с Землей требует значительных вычислительных мощностей и проводить ее для всех астероидов групп Аполлона, Амура и Атона нецелесообразно. После первичного отбора потенциально опасных астероидов для оставшихся астероидов также возможно проведение оценки параметра MOID с начальными данными, вычисленными в результате интегрирования уравнений движения астероидов. Такая процедура дополнительной оценки может применяться для того, чтобы найти астероиды, переходящие в класс потенциально опасных в результате эволюции орбиты. Единственным обязательным требованием метода, предъявляемым к исследуемым телам, является требование конфокальности эллиптических орбит. Введем обозначения для рассматриваемых небесных тел. Пусть имеются небесные тела E и A. Положения небесных тел в пространстве задаются наборами орбитальных элементов: a, e, i, ω, Ω, M , где a - большая полуось, e - эксцентриситет, i - наклонение, ω - аргумент перигелия, Ω - долгота восходящего узла, M - средняя аномалия. Элементы орбит однозначно определяют положение небесного тела. Большая полуось и эксцентриситет задают форму орбиты, наклонение, аргумент перигелия и долготу восходящего узла - ориентацию по отношению к базовой системе координат, средняя аномалия задает положение тела на орбите. В общем случае орбиты тел A и E имеют произвольную ориентацию, то есть параметры i, ω, Ω, M отличны от нуля. Следует отметить, что метод, описанный в данной работе, предназначен для оценки минимального расстояния между орбитами, т. е. данные о реаль146 Быстрая оценка минимального расстояния между . . . ных положениях небесных тел на их орбитах не принимаются во внимание. То есть информация о реальном положении тела на орбите, содержащаяся в орбитальных элементах, в параметре M не используется. Поэтому в дальнейшем вместо шести элементов орбит будет рассматриваться пять. Таким образом, пять элементов орбит будут представлять орбиту небесного тела в пространстве. Чтобы упростить расчеты, осуществим поворот системы координат таким образом, что орбита тела E не будет наклонена. Элементы орбит тела A необходимо будет пересчитать. Отметим, что такое преобразование не изменяет расстояний между орбитами [12]. Преобразование осуществляется с помощью матрицы перехода C, элементы которой зависят от iE , ωE , ΩE : C= cos ΩE cos ωE - sin ΩE cos iE sin ωE cos ΩE sin ωE - sin ΩE cos iE cos ωE sin iE sin ωE sin ΩE cos ωE + cos ΩE cos iE sin ωE sin ΩE sin ωE + cos ΩE cos iE cos ωE - sin iE cos ΩE sin iE sin ωE sin iE cos ωE cos iE . Для расчета новых значений iA , ωA , ΩA необходимо вычислить следующие величины:     x1n x1 c11 + x2 c12 + x3 c13 Xn = x2n  = CX = x1 c21 + x2 c22 + x3 c23  ; x3n x1 c31 + x2 c32 + x3 c33     y1n y1 c11 + y2 c12 + y3 c13 Yn = y2n  = CY = y1 c21 + y2 c22 + y3 c23  ; y3n y1 c31 + y2 c32 + y3 c33     z1n z1 c11 + z2 c12 + z3 c13 Zn = z2n  = CZ = z1 c21 + z2 c22 + z3 c23  , z3n z1 c31 + z2 c32 + z3 c33 где C - матрица перехода, cij - элементы матрицы C; X, Y и Z - векторстолбцы, имеющие следующий вид:   cos ΩA cos ωA - sin ΩA cos iA sin ωA X = sin ΩA cos ωA + cos ΩA cos iA sin ωA  , sin iA sin ωA   - cos ΩA sin ωA - sin ΩA cos iA cos ωA Y = - sin ΩA sin ωA + cos ΩA cos iA cos ωA  , sin iA cos ωA  sin iA sin ΩA Z = - sin iA cos ΩA  . cos iA  Новые значения элементов орбиты тела A рассчитываются по формулам iA = arctan 2 2 z1n + z2n z1n x3n ; ΩA = - arctan ; ωA = - arctan . z3n -z2n y3n Таким образом, перед началом работы алгоритма получим следующий набор орбитальных элементов: для орбиты тела E (aE , eE , 0, 0, 0); для орбиты тела A (aA , eA , iA , ωA , ΩA ), где iA , ωA , ΩA - пересчитанные по приведенным 147 Д е р е в я н к а А. Е. выше формулам элементы орбиты тела A с учетом поворота системы координат [12]. Рассмотрим расположение тел E и A на рис. 1. Солнце - в центре, полярная ось направлена перпендикулярно плоскости орбиты E. Плоскость P содержит полярную ось и проходит через тело A. Таким образом, плоскость P будет перпендикулярна плоскости орбиты E. Так как положение плоскости P зависит от положения тела A на орбите, при движении тела A положение плоскости P будет меняться. Суть описываемого алгоритма заключается в следующем: производится расчет одного оборота тела A по своей орбите с шагом ∆. В результате получается набор расстояний Dv , которые можно рассматривать как дискретную функцию. Затем исходя из имеющихся значений функции устанавливаются области локальных минимумов и производится процедура уточнения. При работе первой части алгоритма положение тела E является зависимым от положения тела A. То есть положение тела E определяется по положению тела A, задаваемому значением истинной аномалии v0 . Положение тела A на орбите задается с помощью истинной аномалии v0 . При этом декартовы координаты тела A и гелиоцентрическое расстояние вычисляются по формулам   a (1 - e2 ) A  rA0 = A ,   1 + eA cos v0  xA0 = rA0 (cos ΩA0 cos(ωA0 + v0 ) - sin ΩA0 sin(ωA0 + v0 ) cos iA0 ), (1)  yA0 = rA0 (sin ΩA0 cos(ωA0 + v0 ) - cos ΩA0 sin(ωA0 + v0 ) cos iA0 ),    z = r (sin(ω + v ) sin i ). 0 A0 A0 A0 A0 Исходя из описанных выше предположений положение тела E на орбите Рис. 1. Расположение орбит тел E и A [Figure 1. Illustration of the orbital geometry of objects E and A] 148 Быстрая оценка минимального расстояния между . . . будет определяться положением тела A. Плоскость P пересекает орбиту тела E в двух точках. Положим, что тело E находится в той точке, что ближе к телу A (см. рис. 1). Таким образом, Солнце, тело E и тело A образуют треугольник, лежащий в плоскости P . Расстояние между E и A обозначим D0 . Найдем координаты тела E, соответствующие положению тела A. Учтем, что в силу вышеупомянутых предположений оба тела находятся вместе с Солнцем в одной плоскости Pv0 . Обозначим долготу тела E как L0 , тогда cos L0 = xA0 , 2 x2 + yA0 A0 sin L0 = yA0 . 2 x2 + yA0 A0 Так как из-за выбора положения орбиты тела E его долгота совпадает с истинной аномалией, гелиоцентрическое расстояние для тела E вычисляется следующим образом: aE (1 - e2 ) E rA0 = . 1 + eE cos L0 Теперь, используя геометрические соображения относительно положений тел, найдем расстояние между телами A и E в их текущих положениях. На рис. 1 этим положениям соответствует расстояние D0 . Рассмотрим треугольник, образованный телами и Солнцем в плоскости Pv0 (см. рис. 2). В треугольнике SAE две стороны равны гелиоцентрическим расстояниям rA0 и rE0 , а третья - искомому расстоянию D0 . Опустим перпендикуляр из точки A на продолжение луча SE, получим точку A . Учитывая, что точка A имеет координаты (xA0 , yA0 , zA0 ), рассчитанные по формулам (1), получим значение SA : SA = 2 x2 + yA0 . A0 Рассмотрим треугольник EA A, в котором EA = SA - SE. Очевидно, EA = 2 (x2 + yA0 )2 - rE0 . A0 Сторона AA равна zA0 . В свою очередь, сторона EA является гипотенузой прямоугольного треугольника, а значит EA2 = EA 2 + AA 2 . Рис. 2. Расположение тел E и A при определении расстояния между орбитами [Figure 2. The triangle SAE for calculating the distance between orbits of E and A] 149 Д е р е в я н к а А. Е. Получается, что искомое значение D0 , соответствующее длине стороны EA треугольника EA A, может быть вычислено следующим образом: 2 2 D0 = zA0 + 2 x2 + yA0 - rE0 A0 2 . После вычисления D0 расчет для одного положения тела A на орбите окончен. Теперь необходимо изменить положение тела A на орбите, увеличив истинную аномалию v0 на величину ∆, и повторить вычисления. Таким образом, получим очередное значение расстояния, например, Dv (см. рис. 1). Проводя описанные выше вычисления до полного оборота тела A по орбите, получим набор расстояний между орбитами тел A и E как дискретную функцию, зависящую от v, для которой можно установить точки локальных минимумов. Один из примеров локального минимума Dmin изображен на рис. 1. Таким образом, первая стадия алгоритма заключается в поиске областей, содержащих локальные минимумы расстояний между орбитами тел E и A. Вторая стадия алгоритма состоит в уточнении локальных минимумов расстояний между орбитами и поиске MOID. Исходя из данных о точках локального минимума и соответствующих им координат тел E и A, зафиксированных на первой стадии работы алгоритма, производится уточнение данных и поиск среди уточненных данных глобального минимума, который и будет являться оценкой значения MOID. Допустим, в результате первой фазы работы алгоритма найдено N предварительных локальных минимумов. Имеем N положений тел E и A и 2N наборов орбитальных элементов. Требуется произвести уточнение позиций E и A для каждого найденного локального минимума Dv . Пусть имеем положения E1 и A1 , а также соответствующее им расстояние между орбитами Dmin 1 . Добавим к рассмотрению еще 4 точки: E1+ , E1- , A1+ , A1- . Эти положения получаются путем варьирования истинной аномалии тела на величину ∆v, т.е. сдвига по направлению движения (индекс +) и против движения (индекс -) на фиксированную величину ∆v. Таким образом, получим 6 позиций и 9 расстояний между орбитами E и A. Расстояния между телами E и A вычисляются как расстояния между двумя точками в декартовой трехмерной системе координат с использованием тривиальных формул:  2  r = aE (1 - eE ) ,  A   1 + eE cos LE  xE = rE cos LE ,  yE = rE sin LE ,    z = 0, E   a (1 - e2 ) A  rA = A ;   1 + eA cos vA  xA = rA (cos ΩA cos(ωA + vA ) - sin ΩA sin(ωA + vA ) cos iA );  yA = rA (sin ΩA cos(ωA + vA ) - cos ΩA sin(ωA + vA ) cos iA );    z = r (sin(ω + v ) sin i ), A A A A A 150 Быстрая оценка минимального расстояния между . . . D2 = (xA - xE )2 + (yA - yE )2 + (zA - zE )2 , где LE и vA - истинные аномалии тел E и A соответственно (с учетом изменений их на величину ∆v). Важно заметить, что на второй итерации цикла количество расстояний, которые необходимо будет вычислить, сократится. К примеру, если новому положению тел E и A будут соответствовать точки E1- , A1- , то добавятся только 2 новых положения тел (т.к. остальные были вычислены на предыдущем шаге). Вторая фаза алгоритма повторяется с учетом того, что положения E1 и A1 равны положениям, соответствующим новому минимуму расстоянию Dmin 1 . Цикл продолжается, пока на k-том шаге цикла найденное минимальное расстояние не совпадет с заданной точностью с найденным на шаге k - 1. Описанный выше алгоритм уточнения выполняется для каждого из N локальных минимумов, найденных на первой стадии работы алгоритма. Вторая стадия повторяется несколько раз с уменьшением шага ∆v. После завершения второй стадии алгоритма минимальное расстояние между орбитами MOID = min Dmin i . i=1,...N Были проведены сравнительные испытания алгоритма с классическим методом G. F. Gronchi, исходный код которого находится в общем доступе (http://adams.dm.unipi.it/?gronchi/kepdist). Начальные данные астероидов брались с сайта Центра малых планет Международного астрономического союза (http://www.minorplanetcenter.net/iau/Dangerous.html), а также с научно-информационного сайта SmallBodies.ru (http://smallbodies. ru). Результаты испытаний приведены в табл. 1 и табл. 2. В табл. 1 приведены результаты оценки количества потенциально опасных астероидов (PHA - Potentially Hazardous Asteroids) среди групп Аполлона, Амура и Атона с помощью метода, описанного в настоящей статье. Астероиды этих групп выбраны в силу того, что почти все потенциально опасные астероиды, известные на сегодняшний день, относятся к одной из этих трех групп. Результаты работы метода совпали с данными, публикуемыми Лабораторией реактивного движения NASA (http://neo.jpl.nasa.gov/orbits) и Центром малых планет международного астрономического союза. Как видно из полученных данных, наибольшее число потенциально опасных астероидов находится в группе Аполлона. Отметим, что при составлении табл. 1 переход Таблица 1 Потенциально опасные астероиды в группах Аполлона, Амура и Атона [Potentially hazardous asteroids found in Apollo, Amor, and Aten groups] Группа астероидов [Group of asteroids] Количество астероидов [Number of asteroids] Количество PHA [Number of PHA] Процент PHA [Percent of PHA] Аполлоны [Apollo] Амуры [Amor] Атоны [Aten] По трем группам [Into three groups] 6 580 4 726 908 1 319 83 145 20.05 % 1.76 % 15.97 % 12 214 1 547 12.67 % 151 Д е р е в я н к а А. Е. Таблица 2 Сравнение точности расчетов MOID [Comparison of the MOID calculation accuracy] Астероид [Asteroid] MOID, а.u. ∆, 10-13 a.u. 2000 SG344 410777 (2009 FD) 99942 Apophis 6344 P-L 2005 WG57 0.000790777 0.002290780 0.000659446 0.028124800 0.001624190 1.5 5.9 10.4 0.7 4.8 астероидов из одной группы в другую во время эволюции орбиты (к примеру, в результате тесных сближений) не рассматривался. Следует отметить, что при сравнении скорости работы в методе, описанном в данной работе, не использовались механизмы параллельных вычислений, которые могли бы еще больше ускорить процедуру расчета MOID. Скорость расчетов составляла в среднем 12 ms на 100 вычислений MOID для описанного алгоритма против 40 ms в среднем на 100 вычислений MOID по методу G. F. Gronchi. В табл. 2 приводятся значения параметра MOID, рассчитанные для нескольких потенциально опасных астероидов. В силу того, что при расчете параметра MOID по методу, описанному в данной статье, получаются значения, совпадающие с классическим методом G. F. Gronchi как минимум до 9 разряда, в табл. 2 во втором столбце приведено значение MOID с точностью до 9-го знака, а в третьем - разность ∆ = |MOIDGr - MOIDFast |, где MOIDGr - значение, полученное по методу G. F. Gronchi, a MOIDFast - значение, полученное по описанному в статье алгоритму. Таким образом, видно, что в результате настройки количества итераций в первой и второй частях алгоритма (параметры ∆ и ∆v) можно добиться совпадения значений MOID с рассчитанными по методу G. F. Gronchi до значений порядка от 10-12 до 10-14 . Характеристики рабочей станции, на которой проводились сравнительные испытания методов: CPU: Intel Core i7-4800MQ (4 kernels, 2.7 GGz); RAM: 8.0 GiB, DDR3; HDD: 1 TiB, 7200 rpm.; OS: Windows 8.1 64 bit. Метод, описанный в данной работе, может применяться не только для определения минимального расстояния между орбитами Земли и астероида, но и для расчета минимального расстояния между орбитами двух произвольных небесных тел. Как уже было отмечено выше, единственным условием остается требование конфокальности эллиптических орбит двух небесных тел. Представленный метод расчета MOID, в отличие от аналитических методов К. В. Холшевникова и G. F. Gronchi, не предоставляет информацию об особых точках функции, описывающей расстояние между орбитами E и A (максимумы расстояний, седловые точки). Однако так как метод создан для оценки только параметра MOID для большого количества астероидов с высокой скоростью и точностью, данная информация для поставленных задач не требуется.

About the authors

Andrey E Derevyanka

Samara State Technical University

Email: AndrDerev@gmail.com
Postgraduate Student, Dept. of Applied Mathematics & Computer Science 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

References

Supplementary files