Астрономический вестник, 2020, T. 54, № 6, стр. 567-576

ВВЕДЕНИЕ

В 2004 г. в обсерватории Китт‒Пик в Аризоне был открыт астероид 2004 MN4, впоследствии получивший имя Апофис и номер 99942. Крупный размер, около 300 м, и существенная вероятность соударения с Землей вызвали интерес многих исследователей, в том числе в СПбГУ. Особенностью данного астероида является тесное сближение с Землей в 2029 г., порождающее множество резонансных возвратов и затрудняющее точное предсказание дальнейшей траектории (Соколов и др., 2008). Хотя к настоящему моменту соударения вскоре после 2029 г. исключены, сохраняется вероятность соударения во второй половине XXI века. https://cneos.jpl.nasa.gov/sentry/ – результаты расчета вероятности соударения, полученные NASA.

Для поиска и изучения возможных соударений на кафедре небесной механики СПбГУ был разработан программный комплекс v19. С его помощью для ряда астероидов был получен существенно более полный список возможных соударений, чем приводится на сайте NASA, в том числе для астероида Апофис (Соколов и др., 2013). Поиск соударений осуществляется путем одномерного варьирования начальных данных. Типичное время работы над одним астероидом составляет от нескольких часов до нескольких суток. Расширение исследований на большее количество астероидов и переход к многомерному варьированию потребовали бы колоссальных затрат вычислительных ресурсов.

Основное время занимает вычисление траектории астероида и параметров сближений с Землей. Для этой задачи было решено создать с нуля новую программу. Программа R⁰ предназначена для расчета траекторий и сближений с планетами, Луной и Солнцем большого числа виртуальных астероидов. Модель движения учитывает притяжение больших планет, координаты которых берутся из эфемерид DE430. Для интегрирования уравнений выбран метод Гаусса–Эверхарта (Авдюшев, 2010). Создание с нуля позволило провести существенную оптимизацию. Ускорение по сравнению с программным комплексом v19 составило более 10 раз. Точное число зависит от конкретного теста. Новая программа опробована на задаче определения вероятности соударения методом Монте‒Карло. Полученные результаты близки к приведенным на сайте NASA.

УРАВНЕНИЕ ДВИЖЕНИЯ АСТЕРОИДА

Пусть $x$ – вектор положения астероида, ${{x}_{i}}$ – вектор положения массивного тела Солнечной системы, ${{r}_{i}}$ – расстояние между массивным телом и астероидом, $G$ – гравитационная постоянная, ${{m}_{i}}$ – масса тела. Тогда уравнение движения можно записать в следующем виде:

(1)

Слагаемое $F$ – позволяющее учесть различные возмущения, принято равным нулю в рамках данной работы. Преобразовав уравнение движения, для трехмерного пространства можно получить систему из 6 обыкновенных дифференциальных уравнений первого порядка. Начальными данными будут координаты и скорости в стартовый момент времени. Для реального астероида начальные данные определены с погрешностями, так что существует область неопределенности ненулевого размера, где может быть астероид. Интерес представляют области начальных данных, точки из которых ведут к соударению с Землей или другими телами. Пересечение таких областей с областью неопределенности означает ненулевую вероятность соударения.

МЕТОД ГАУССА–ЭВЕРХАРТА

Подробное описание метода можно найти в работе (Авдюшев, 2010), здесь же будут изложены основные формулы и алгоритм интегрирования на шаге. Рассмотрим следующую задачу:

Здесь $t$ – независимая переменная, $x$ – интегрируемая переменная, $f$ – заданная функция. На шаге интегрирования $h$ введем переменную $\tau = {{\left( {t - {{t}_{0}}} \right)} \mathord{\left/ {\vphantom {{\left( {t - {{t}_{0}}} \right)} h}} \right. \kern-0em} h}$, а правую часть аппроксимируем многочленом степени $k$:

Интегрируя многочлен, получаем:

Перепишем правую часть в виде интерполяционного многочлена Ньютона на сетке ${{\tau }_{0}},{{\tau }_{1}}, \ldots ,{{\tau }_{k}},{{\tau }_{0}} = 0$:

Если известны коэффициенты $A$ на предыдущем шаге, начальным приближением на новом шаге будет:

r – отношение нового шага интегрирования к предыдущему; чаще всего r = 1.

Коэффициенты $e$ вычисляются по рекуррентным формулам:

Прибавив полученную на предыдущем шаге разницу $\Delta A$ между $A$ в конце итерационного процесса и начальным приближением $\bar {A}$, оценку можно улучшить. В начале интегрирования за начальное приближение принимаются нулевые значения.

Итерация производится следующим образом. Для каждого $i = 1,2, \ldots ,k$ последовательно выполняется цикл $A \to {{x}_{i}} \to {{f}_{i}} \to {{\alpha }_{i}} \to A$. Первое действие $A \to {{x}_{i}}$ производится по формуле (4); второе действие ${{x}_{i}} \to {{f}_{i}}$ – в соответствии с уравнением (1). Для третьего действия ${{f}_{i}} \to {{\alpha }_{i}}$ преобразуем формулу (5):

За одно прохождение цикла вычисляется одно значение ${{\alpha }_{i}}$. Начальное приближение $\alpha $ выводится из начального приближения $A$:

Четвертое действие $\alpha \to A$ будет обратным к (9):

Коэффициенты $c$ и $d$ вычисляются по рекуррентным формулам:

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

Для повышения порядка метода узлы ${{\tau }_{i}}$ следует выбирать из разбиений Гаусса–Лежандра, Гаусса–Радо или Гаусса–Лобатто. В программе R⁰ использовано разбиение Гаусса–Лобатто, узловые значения для которого являются корнями уравнения:

Порядок метода в этом случае равен $2k$, ${{\tau }_{0}} = 0,{{\tau }_{k}} = 1$.

ОПТИМИЗАЦИЯ МЕТОДА ДЛЯ КОНКРЕТНОЙ ЗАДАЧИ

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

Оценим количество арифметических операций, затрачиваемое на каждое из четырех действий в течение одной итерации при $k = 6$. Для одного астероида интегрируемая переменная представляет собой 6-вектор, но, например, при вычислении правых частей существенны только 3 компоненты из 6, а расстояние до планеты достаточно вычислить один раз, как и величину ${{G{{m}_{i}}} \mathord{\left/ {\vphantom {{G{{m}_{i}}} {{{r}_{i}}^{3}}}} \right. \kern-0em} {{{r}_{i}}^{3}}}$. Для первого действия $A \to {{x}_{i}}$ величины ${{\tau _{i}^{{j + 1}}} \mathord{\left/ {\vphantom {{\tau _{i}^{{j + 1}}} {\left( {j + 1} \right)}}} \right. \kern-0em} {\left( {j + 1} \right)}}$ лучше вычислить заранее, тогда будет затрачено $6k\left( {k + 1} \right) = 252$ сложения и $6k\left( {k + 1} \right) = 252$ умножения. Второе действие ${{x}_{i}} \to {{f}_{i}}$ для 10 массивных тел потребует $30k = 180$ вычитаний, $10k = 60$ делений, $10k = 60$ извлечений квадратного корня, $30k = 180$ сложений и $80k = 480$ умножений. Для третьего действия ${{f}_{i}} \to {{\alpha }_{i}}$, поскольку операция деления дороже умножения, имеет смысл вычислить величины ${1 \mathord{\left/ {\vphantom {1 {\left( {{{\tau }_{i}} - {{\tau }_{j}}} \right)}}} \right. \kern-0em} {\left( {{{\tau }_{i}} - {{\tau }_{j}}} \right)}}$. Тогда третье действие потребует всего $6k{{\left( {k + 1} \right)} \mathord{\left/ {\vphantom {{\left( {k + 1} \right)} 2}} \right. \kern-0em} 2} = 126$ вычитаний и $6k{{\left( {k + 1} \right)} \mathord{\left/ {\vphantom {{\left( {k + 1} \right)} 2}} \right. \kern-0em} 2} = 126$ умножений. Четвертое действие $\alpha \to A$ выполняется за $6{{k}^{2}}{{\left( {k + 1} \right)} \mathord{\left/ {\vphantom {{\left( {k + 1} \right)} 2}} \right. \kern-0em} 2} = 756$ умножений и $6{{k}^{2}}{{\left( {k - 1} \right)} \mathord{\left/ {\vphantom {{\left( {k - 1} \right)} 2}} \right. \kern-0em} 2} = 540$ сложений. Последнее действие, как можно заметить, сопоставимо по затратам с вычислением правых частей. При учете дополнительных сил соотношение изменится в пользу правых частей, но при повышении порядка разбиения $k$ затраты на четвертое действие растут очень быстро. Например, при $k = 8$ потребуется уже 1728 умножений и 1344 сложения на одну итерацию. Можно было бы оптимизировать это действие, но можно объединить с первым в сокращенной схеме $\alpha \to {{x}_{i}} \to {{f}_{i}} \to {{\alpha }_{i}}$. Тогда получение начального приближения будет осуществляться по схеме $\alpha \to A \to \bar {A} \to \bar {\alpha }$, а вместо поправки $\Delta A$ можно использовать аналогичную $\Delta \alpha $. Объединенное действие $\alpha \to {{x}_{i}}$ выполняется по формуле:

Количество арифметических операций такое же, как и для действия $A \to {{x}_{i}}$, при этом четвертое действие полностью отсутствует. Аналогичным образом упрощается получение начального приближения:

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

Наконец, поскольку исходная система уравнений имеет второй порядок, можно использовать другой вариант аппроксимации решения (Everhart, 1974) и таким образом избавиться от необходимости вычислять скорости в течение итерации:

(18)

Тогда действие $\alpha \to {{x}_{i}}$ примет вид:

Основной идеей оптимизации является параллельное интегрирование большого числа астероидов, поэтому информация о них всех занимает много места в памяти. При тестировании оказалось, что поправка $\Delta \alpha $ не приводит к существенному росту производительности, зато повышает расход памяти почти в два раза. Начальное приближение $\bar {\alpha }$ уже без поправки позволяет делать всего 2–3 итерации на шаге. В связи с этим программа R⁰ эту поправку не использует.

Итерации производятся до полной сходимости положения и скорости на конце шага. Критерий выбран с учетом машинной точности, составляющей 52 двоичных порядка для типа double.

АЛГОРИТМ ВЫБОРА ШАГА И ГРУППИРОВКИ ВИРТУАЛЬНЫХ АСТЕРОИДОВ

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

В качестве базового шага выбрана величина ${{h}_{0}} = {{4}^{d}}$, обеспечивающая наиболее удобное использование эфемерид DE430. В течение базового шага сохраняется аппроксимация положений планет, Луны и Солнца. Ни один шаг интегрирования не превышает этой величины, а большинство равны ей. Критерий уменьшения шага – эвристический, направленный прежде всего на обнаружение сближений. Для уменьшения шага хотя бы одно отношение планетоцентрической скорости к планетоцентрическому расстоянию должно вырасти до некоторой регулируемой величины $2M$. Если отношение превышает ${{2}^{i}}M$, шаг уменьшается до ${{2}^{{ - i}}}{{h}_{0}}$. В итоге шаг подбирается из расчета, чтобы изменение планетоцентрического радиус-вектора не превысило четверти величины этого радиус-вектора (или 1/25 длины окружности). Порядок разбиения $k$ подобран для достижения наилучшей точности и равен 5. Точность контролировалась параллельным интегрированием с меньшим шагом.

В ходе интегрирования на базовом шаге производятся следующие итерации: определяется множество астероидов с наименьшим значением времени, до которого рассчитана траектория; для каждого из них определяется шаг; затем для каждого шага выполняется расчет эфемерид и собственно интегрирование на шаге. Если шаг можно увеличить по сравнению с предыдущим, проверяется текущая степень дробления шага. Например, если астероид был рассчитан на 5/8 базового шага, новый шаг выбирается не более 1/8, чтобы в следующий раз он смог войти в группу, интегрируемую с шагом 1/4.

АЛГОРИТМ ПОИСКА И УТОЧНЕНИЯ СБЛИЖЕНИЙ

Для каждой планеты задается расстояние ${{R}_{{{\text{min}}}}}$, при сближении до которого нужно вычислить параметры сближения. После каждого шага интегрирования проверяется наличие сближения в пределах шага. Для этого вычисляется скалярное произведение планетоцентрических скорости и положения астероида в начале и в конце шага. Если в начале произведение отрицательно, а в конце положительно, произошло сближение. Если расстояние на конце шага меньше $2{{R}_{{{\text{min}}}}}$, производится уточнение сближения. Если полученное минимальное расстояние оказывается меньше ${{R}_{{{\text{min}}}}}$, параметры сближения записываются в файл: номер виртуального астероида, дата сближения, номер тела, минимальное расстояние. Если минимальное расстояние меньше радиуса тела, то сближение является соударением, расчет траектории в этом случае останавливается.

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

СРАВНЕНИЕ РЕЗУЛЬТАТОВ R⁰ И v19

Первое требование к программе – точность результата. Как бы быстро программа не считала, если результат отличается от контрольного больше, чем на радиус Земли, для предсказания соударений такая программа не годится по меньшей мере в условиях данного теста. Для программы R⁰ контрольный результат предоставила старая программа v19, использующая ту же модель движения, но другой интегратор. Основным тестом стал расчет астероида (1036) Ганимед на 320 000 суток. Большое время интегрирования (в 10 раз больше обычного) и наличие тесного сближения (с Марсом) делали задачу сложнее или сопоставимой с предположительными задачами для большинства АСЗ. Расхождение положения астероида в течение всего времени интегрирования оставалось менее 10^–6 а. е.

Второе требование к программе – скорость работы. Вообще говоря, старая программа не может одновременно рассчитывать сближения с разными планетами и должна быть запущена отдельно для каждой из планет. Уже поэтому она заведомо в 10 раз медленнее новой. Однако в задаче поиска сближений и соударений конкретно с Землей этот недостаток не имеет значения. Для теста была взята тысяча виртуальных астероидов в окрестности номинальной орбиты Апофиса и процессор Intel Pentium G3260. Расчет велся на 20 000 суток вперед (около 55 лет). Среднее время, затраченное на один астероид, составило у программы v19 – 300 мс, у программы R⁰ – 13 мс.

ВЫЧИСЛЕНИЕ ВЕРОЯТНОСТИ СОУДАРЕНИЯ МЕТОДОМ МОНТЕ‒КАРЛО

Данная задача требует расчета траекторий большого числа виртуальных астероидов, что позволит максимально использовать преимущества новой программы. С помощью программы R⁰ был произведен расчет траектории до 2132 года 10⁷ виртуальных астероидов из области неопределенности, выбранных случайно в соответствии с нормальным распределением, задаваемым математическим ожиданием и матрицей ковариаций элементов орбиты. https://ssd.jpl.nasa.gov/sbdb.cgi – база данных малых тел Солнечной системы. Для генерации случайных чисел использован вихрь Мерсенна. Астероиды выбраны из числа имеющих размер более 30 м и вероятность соударения с Землей более 10^–8 по данным NASA, с которыми и сравниваются результаты.

В табл. 1 приведены вероятность соударения с Землей по данным NASA и количество попавших в тела Солнечной системы виртуальных астероидов из 10⁷ по расчетам программой R⁰. Для удобства сравнения вероятность NASA умножена на 10⁷. Всего было исследовано 200 астероидов.

ЗАКЛЮЧЕНИЕ

Создана программа R⁰, значительно ускорившая расчет траекторий опасных астероидов. Различие результатов со старой программой v19, а также с результатами NASA, невелико, что показывает потенциальную пригодность программы R⁰ в задачах астероидно-кометной безопасности. Хотя, например, астероиды 29 075 и 101 955 (в конце таблицы) демонстрируют существенное различие, оно естественным образом проистекает из отсутствия возможных соударений до 2132 г. и их наличием значительно позже. В большинстве же случаев интервал исследования NASA меньше, поэтому вероятность по результатам R⁰ обычно больше. Следует, однако, с осторожностью полагаться на полученные вероятности. В ряде случаев положенные в основу предположения и ошибки самого метода могут исказить результат на порядок и более.

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

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

Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 19-32-90149.