7.1. Определение предварительной орбиты и ее последующие уточнения. Оценка точности элементов орбиты
Для выделения потенциально опасных астероидов из общего числа АСЗ, для оценки вероятности столкновения их с Землей и предотвращения столкновений первостепенное значение имеют знание параметров движения и оценка их вероятных ошибок.
Как известно, движение тела относительно некоторой инерциальной системы координат полностью определяется действующими на него силами и начальными условиями. В качестве последних обычно выбирают координаты и компоненты скорости в некоторый момент времени или шесть элементов орбиты. Обратная задача заключается в том, чтобы по наблюдаемому движению небесного тела определить начальные условия движения, например элементы орбиты в некоторый момент времени. Так как каждое позиционное наблюдение дает две сферические координаты (α и δ), то минимальное количество наблюдений, необходимых для определения шести элементов эллиптической орбиты, равно трем. Орбита, найденная по трем или небольшому числу наблюдений, называется предварительной. Для определения предварительной орбиты большей частью используются методы, основанные на работах Лагранжа, Гаусса и Лапласа [Субботин, 1968; Херрик, 1977; Быков, 1989; Marsden, 1991]. Как правило, при определении предварительной орбиты астероида или кометы учитывается только притяжение Солнца. Возмущающим влиянием больших планет и другими возможными возмущениями в движении тела при этом пренебрегают.
Предварительная орбита имеет невысокую точность как из-за ошибок наблюдений, на основе которых она определена, так и из-за пренебрежения действующими на тело силами. Однако определение предварительной орбиты является необходимым этапом, поскольку оно позволяет вычислить эфемериду тела для продолжения наблюдений в ближайшие дни и не потерять объект. Если в дальнейшем удается провести дополнительные наблюдения или найти в каталогах наблюдения, принадлежащие тому же телу, то предварительная орбита подвергается исправлению, или уточнению, с учетом старых и новых наблюдений. При этом уже учитываются возмущения, вызываемые другими телами Солнечной системы помимо Солнца, и, возможно, иные возмущения.
Уточнение параметров движения чаще всего выполняется по методу наименьших квадратов (МНК). Напомним основные положения этого метода и некоторые формулы, используемые в дальнейшем.
Вновь наблюдаемые координаты тела, как правило, заметно отличаются от тех координат, которые вычисляются согласно теории движения с первоначально найденными параметрами (элементами орбиты тела). Процесс уточнения предварительной орбиты сводится к тому, чтобы найти такие поправки к исходной системе элементов, которые уменьшали бы рассогласование между наблюденными и вычисленными положениями.
Пусть имеются n наблюденных положений тела, которые обозначим как Ok (под Ok можно понимать наблюденное значение любой координаты). По теории движения с исходной системой элементов орбиты E0i (i = 1…, 6) на моменты наблюдений tk вычисляются положения Ck:
F (t k ; E 0 1 …, E 0 6 ) = Ck.
Разности Ok — Ck (их обычно называют «наблюденное минус вычисленное»), с одной стороны, зависят от неточности исходной системы элементов, а с другой, — от ошибок наблюдений, причем вклад первой составляющей на первых порах оказывается существенно больше второй. Функцию F в окрестности исходного значения F (t; E0i) можно представить по степеням приращений элементов орбиты:
Предположим, что ошибки наблюдений малы, и что E1…, E6 есть та система элементов, которая позволяет более точно вычислить наблюдаемые значения Ok. Если допустить, что она мало отличается от исходной системы E0i, и что высшими степенями приращений ΔE0i можно поэтому пренебречь, то формула (7.1) позволяет написать
Это так называемое условное уравнение относительно искомых поправок ΔE0i.
Частные производные в левых частях условных уравнений могут считаться известными функциями, поскольку они всегда могут быть вычислены с большей или меньшей точностью.
Поскольку число наблюдений при уточнении орбиты почти всегда больше числа уточняемых параметров, то система условных уравнений является избыточной. В общем случае речь может идти лишь о ее приближенном решении. В методе наименьших квадратов решение ищется на основе принципа Лежандра — минимизации суммы квадратов остаточных уклонений. Под остаточными уклонениями εk понимаются разности между левыми и правыми частями уравнений (7.2):
Согласно принципу Лежандра, искомые неизвестные поправки должны минимизировать величину
Необходимые условия минимума S как функции переменных ΔE0i записываются в виде
Эти условия образуют систему из шести линейных уравнений относительно шести неизвестных ΔE0i (i = 1…, 6). Например, первое из них записывается в виде
Остальные уравнения записываются аналогично.
Система из шести уравнений (7.3) относительно неизвестных ΔE0i называется нормальной. Использование матричного исчисления позволяет представить нормальную систему и ее решение в компактном виде. Составим матрицу коэффициентов условных уравнений:
Обозначим также вектор-столбец с компонентами ΔE0i как вектор X, а вектор-столбец правых частей с компонентами Ok — Ck как вектор L. В таком случае система условных уравнений запишется в виде
BX = L.
Нормальная система записывается в виде
B T BX = B T L,
где символ T означает транспонирование матрицы (замену строк столбцами и наоборот).
Симметричную матрицу нормальной системы BTB обозначим буквой Q. Решение нормальной системы может быть найдено умножением обеих частей уравнения на матрицу
Q -1 = (B T B) -1 ,
где символом Q-1 обозначена матрица, обратная матрице Q (заметим, что матрица Q-1 как обратная симметричной матрице, также является симметричной). Произведение обратной матрицы на саму матрицу дает единичную матрицу, вследствие чего решение записывается в виде
X = Q -1 B T L. (7.4)
Складывая найденные поправки ΔE0i с исходной системой параметров E0i, находят исправленную систему. Поскольку при образовании системы условных уравнений мы пренебрегли высшими степенями поправок, то исправленная система элементов не обеспечивает минимального значения суммы квадратов остающихся невязок, хотя обычно уменьшает ее. Для достижения минимума процедуру дифференциального исправления системы элементов приходится повторять до тех пор, пока поправки к элементам не станут достаточно малыми. Найденное таким образом решение называют номинальным.
На практике используется большое число методов решения нормальной системы, в том числе и тот, который, согласно (7.4), основан на обращении матрицы Q, хотя его следует избегать в случае малости определителя матрицы. В теоретическом плане представление решения в виде (7.4) является наглядным и позволяет раскрыть ряд особенностей этого решения. К этому вопросу мы еще вернемся, но прежде рассмотрим вероятностный смысл решения системы условных уравнений методом МНК.
Интересующие нас особенности имеют место только в случае нормального закона распределения ошибок (закона Гаусса):
(e — основание натурального логарифма).
Центр распределения определяется значением элемента в номинальном решении, а дисперсии элементов определяются диагональными элементами матрицы Q-1 — обратной матрицы нормальной системы. Иначе говоря, если обозначить среднеквадратичную ошибку элемента Ei как σi, то
где σ — средняя квадратичная величина остаточных уклонений:
n — число условных уравнений, m — число определяемых неизвестных, в нашем случае 6, Q-1ii — i-й диагональный элемент обратной матрицы нормальной системы.
Для нормального распределения ошибок элементов орбиты справедливы особенности нормального распределения, в частности то, что вероятность появления ошибки, превышающей утроенное значение среднеквадратичной ошибки элемента, меньше 0,003.
Когда мы говорим об ошибках элементов орбиты, то понимаем при этом возможные отличия элементов от тех значений, которые они имеют в номинальном решении. Таким образом, областью возможных значений для каждого элемента является Ei ± 3σi. Каждую возможную орбиту можно представить как точку шестимерного пространства, по осям которого откладываются значения элементов орбит. Рассмотрим малую окрестность некоторой точки этого пространства. Вероятность попадания орбиты в эту окрестность зависит от одновременного попадания шести элементов орбиты в соответствующие элементарные интервалы ΔEi. Мы уже видели, что при нормальном распределении ошибок эти вероятности определяются формулами типа (7.5), т. е.
где под xi следует понимать элемент Ei, p(xi) — плотность вероятности распределения ошибок соответствующего элемента, σi — корень квадратный из дисперсии ошибок (среднеквадратичная ошибка i-го элемента).
Предположим ради простоты изложения, что случайные ошибки элементов Ei и Ej попарно независимы, т. е. вероятность попадания ошибки элемента Ei в некоторый интервал не зависит от ошибки элемента Ej. В этом случае ошибки всех элементов являются независимыми в совокупности. Плотность вероятности одновременного попадания шести элементов в достаточно малую окрестность точки (E1…, E6) в этом случае выражается как произведение плотностей вероятностей распределения ошибок отдельных элементов:
p(E 1 …, E 6 ) = p 1 (E 1 )p 2 (E 2 )… p 6 (E 6 ).
Каждый сомножитель в правой части последней формулы определяется формулой типа (7.5). Из этого вытекает, что плотность вероятности в точке r в случае шестимерного нормального распределения при сделанном предположении определяется формулой
Указанная плотность вероятности остается неизменной во всех точках пространства, где
При любом положительном значении постоянной это выражение представляет собой уравнение эллипсоида в осях, совпадающих по направлению с главными осями эллипсоида и имеющих начало в точке (x10, x20, x30…, x60) шестимерного пространства. Если представить, что в шестимерном пространстве элементов по осям прямоугольной системы координат с началом в точке, отвечающей номинальной орбите, отложены величины σi и представить себе шестимерный эллипсоид с полуосями σi, то плотность вероятности на таком эллипсоиде будет всюду одинаковой. То же самое будет справедливо и для любого другого подобного и подобным образом расположенного эллипсоида. Такие эллипсоиды называются эллипсоидами равных плотностей вероятностей.
По аналогии с одномерным случаем можно заключить, что вероятность попадания точки внутрь некоторого эллипсоида равна интегралу
где интегрирование распространяется на все пространство, ограниченное эллипсоидом. Если полуоси эллипсоида неограниченно увеличиваются, то интеграл по всему пространству равен единице. Если представить эллипсоид с полуосями, равными 3σi, то вероятность попадания точки в область пространства, ограниченную этим эллипсоидом, близка к единице (0,99736). Такой эллипсоид будем называть доверительным.
Выше предполагалось, что ошибки элементов независимы. На самом деле они корреляционно связаны. Отражением этих связей между ошибками отдельных элементов, найденных по методу МНК, являются величины недиагональных элементов обратной матрицы Q-1, которую называют корреляционной матрицей решения или матрицей ковариаций. Корреляционные связи могут проявляться по-разному. Примером двух элементов, находящихся в жесткой корреляционной зависимости, являются долгота узла и угловое расстояние перигелия от узла при малом наклоне орбиты. Ошибки этих величин близки по величине и противоположны по знаку.
Сделанное выше допущение о независимости случайных ошибок элементов эквивалентно допущению, что все недиагональные элементы матрицы ковариаций равны нулю. В том случае, если это допущение неверно, плотность вероятности многомерного нормального распределения будет иметь более сложный вид по сравнению с (7.7). В показателе экспоненты будет присутствовать сумма не только квадратов, но и смешанных членов вида (xi — xi0)(xj — xj0) с коэффициентами, зависящими от недиагональных элементов матрицы ковариаций (коэффициентов корреляции). Приравнивание суммы в показателе экспоненты к положительной постоянной дает уравнение эллипсоида равной плотности вероятности, но в этом случае ориентация главных осей эллипсоида не совпадает с ориентацией координатных осей. Путем поворота координатных осей уравнение эллипсоида может быть приведено к виду (7.8), в котором отсутствуют смешанные члены.
Корреляционные матрицы, определяющие погрешности элементов и корреляционные связи между ними, находят важное применение при определении погрешностей различных функций этих элементов. Этот вопрос еще будет обсуждаться в следующих параграфах.
Подводя итог, важно обратить внимание на то, что элементы истинной орбиты тела остаются неизвестными. Любая точка внутри доверительного эллипсоида представляет некоторую орбиту, совместимую с имеющимися наблюдениями. Однако вероятность того, что реальная орбита находится в малой окрестности номинального решения, является максимальной по сравнению с другими возможными решениями.
Отметим, что до сих пор мы рассматривали все наблюдения как имеющие одинаковую точность. На практике приходится определять элементы орбиты на основе рядов наблюдений, выполненных с различными точностями (имеющими различные среднеквадратичные ошибки σ1, σ2…, σn). В таких случаях вводят понятие веса наблюдения, определяя его как
где σ0 — произвольное положительное число.
Решение системы условных уравнений в таком случае ищут исходя из обобщенного принципа Лежандра: решение системы должно минимизировать взвешенную сумму квадратов остающихся невязок:
Из этого требования вытекает правило преобразования системы условных уравнений и ее решения: каждое условное уравнение должно быть умножено на корень квадратный из веса соответствующего наблюдения. После этой операции (так называемого приведения к равноточным наблюдениям) система решается так же, как в случае наблюдений, имеющих одну и ту же среднюю ошибку.
7.2. Нелинейный характер распространения ошибок начальных данных. Поиск потенциально опасных сближений астероидов с Землей и оценка вероятности столкновений
После того как номинальная орбита астероида определена, появляется возможность предвычислить его движение в предстоящий период времени и определить, угрожает ли Земле столкновение с ним в обозримом будущем. В зависимости от точности найденной орбиты такие расчеты желательно выполнять для всех АСЗ на интервалах от нескольких лет до нескольких десятков лет, а иногда и до нескольких сотен лет. Прогнозирование движения выполняется методом численного интегрирования уравнений движения, в которых учитываются члены, обусловленные притяжением больших планет и наиболее массивных астероидов (в случаях, требующих особой точности, иногда учитываются возмущения от трехсот наиболее массивных астероидов, см. раздел 7.3). В ходе численного интегрирования фиксируются моменты тесных сближений с Землей и другими большими планетами, которые могут заметным образом трансформировать орбиту тела и тем самым оказать влияние на ее последующие сближения с Землей.
Поскольку столкновения достаточно крупных тел с Землей — весьма редкие события, то при прогнозировании движения тела по номинальной орбите столкновения с Землей, как правило, не обнаруживаются. Нужно, однако, иметь в виду, что номинальная орбита является лишь одной из бесчисленного количества других возможных орбит, элементы которых более или менее близки к элементам номинальной орбиты. Фактическая орбита тела, которая нам не известна, находится где-то внутри области, ограниченной доверительным эллипсоидом (см. раздел 7.1).
Аналогичное представление об области неопределенности начальных условий движения справедливо и в том случае, если рассматривать точки не в пространстве элементов орбит, а в пространстве начальных значений прямоугольных координат и скоростей тел, что имеет несколько большую наглядность.
По мере увеличения числа использованных наблюдений и расширения покрываемого ими временного интервала, ошибки определения элементов, вообще говоря, уменьшаются, а вместе с тем сокращаются и полуоси доверительного эллипсоида. Его центр, соответствующий новому номинальному решению, при этом также несколько смещается в пространстве.
Каждая точка внутри доверительного эллипсоида соответствует некоторой возможной орбите. Тело на возможной орбите мы будем называть виртуальным (возможным) астероидом [Milani et al., 2002].
Если внутри доверительного эллипсоида случайным образом выбрать большое число виртуальных астероидов и следить за их движением на некотором интервале времени, то можно наблюдать, как с течением времени изменяются форма и размеры области пространства, в которой в данный момент заключены виртуальные астероиды. Во всех случаях, с которыми приходится иметь дело на практике, область, первоначально занятая доверительным эллипсоидом, постепенно расширяется и вытягивается вдоль номинальной орбиты тела. Причиной этого являются небольшие различия элементов орбит виртуальных астероидов, причем различие в среднем движении вызывает пропорциональные времени расхождения в средней аномалии, значение которой определяет положение тела на орбите. В результате граница области, занятой виртуальными астероидами, постепенно превращается в очень вытянутый эллипсоид, который можно представить в виде трубки более или менее постоянной ширины, окружающей номинальную орбиту. С течением времени длина трубки может достичь тысяч и миллионов километров и даже превзойти длину орбиты тела.
Большие искажения области пространства, занятой виртуальными астероидами, обусловливают их тесные сближения с Землей или другими планетами. Орбиты с близкими начальными условиями движения по прошествии большого интервала времени могут оказаться весьма далекими друг от друга или, напротив, скрещивающимися друг с другом, что может быть квалифицировано как наложение области возможных движений самой на себя. Во всех этих случаях принято говорить о нарушении линейности задачи. Математически это означает, что приращение некоторой функции начальных значений параметров существенно отличается от ее первого дифференциала и при ее вычислении нельзя пренебрегать членами с дифференциалами высших порядков.
Решение задачи об оценке вероятности встречи астероида с Землей мы рассмотрим, следуя в целом линии, намеченной в работах [Milani et al., 2000; 2002]. На первом этапе будем предполагать, что задача имеет линейный характер, отложив на потом более сложные случаи. Фактически это равносильно предположению, что область пространства, занятая виртуальными астероидами в окрестности сближения номинальной орбиты с Землей, представляет собой эллипсоид, хотя его размеры и форма (вытянутость) отличаются от размера и формы доверительного эллипсоида в начальный момент времени.
7.3. Возмущения, которые необходимо учитывать при уточнении орбиты и прогнозе движения
Вычисление возможности столкновения того или иного небесного тела с Землей или иной планетой путем достаточно точного прослеживания траектории его движения на длительном интервале времени является одной из сложных и трудоемких задач вычислительной астрономии. Задача требует целесообразно полного учета действующих на тело сил и достаточно точного метода интегрирования уравнений его движения, в которых учтены все принимаемые во внимание силы. В совокупности эти факторы составляют модель движения тела. Учитываемые факторы в модели движения должны сообразовываться друг с другом. Бессмысленно учитывать в уравнениях движения малые по величине силы, если метод численного интегрирования не в состоянии обеспечить необходимую точность вычислений на всем интервале прогноза. С другой стороны, совершенно излишне использовать очень точный метод численного интегрирования, если действующие на тело силы неизвестны достаточно точно и результат влияния этих сил может на много порядков превзойти ошибку, зависящую от погрешности численного интегрирования. Прогноз движения должен также сообразовываться со знанием начальных условий движения, которые, как правило, определяются из наблюдений, обремененных теми или иными ошибками. Без учета возможных ошибок результаты прогноза могут оказаться ложными или неполными.
Для большей конкретики рассмотрим этот вопрос на примере учета светового давления в движении астероида (99942) Апофис. Световое давление оказывает заметное влияние на движение этого сравнительно небольшого астероида (D ∼ 270 м [Cellino et al., 2007]). Приближенно учет светового давления может быть выполнен по известным формулам [Аксенов, 1977], если форму поверхности астероида считать сферической и принять для него некоторые значения радиуса, плотности и коэффициента диффузного отражения поверхности.
При включении светового давления в число возмущающих факторов на двухгодичном интервале с 15 марта 2004 г. по 16 августа 2006 г. в ходе уточнения орбиты из наблюдений была найдена система элементов орбиты Апофиса. Прогнозирование движения астероида вперед на 25 лет с данной системой элементов при сохранении модели движения неизменной дает минимальное расстояние Апофиса от Земли 13 апреля 2029 г., равное 38 220 км. Если при выполнении прогноза не учитывать световое давление, то результат оказывается почти на 80 км меньше. Результаты этих вычислений понятны. Позиционные наблюдения астероида требуют учета влияния светового давления. При включении светового давления в число возмущающих факторов большая часть эффекта светового давления на наблюдения учитывается по формулам, даже если форма поверхности, масса и коэффициент диффузного отражения известны не вполне точно. Оставшаяся часть эффекта в некоторой степени учитывается подгонкой параметров орбиты к наблюдениям при их определении по методу МНК. Таким образом, включение светового давления в число возмущающих факторов является необходимым, если желательно обеспечить километровую точность прогноза.
Какие еще возмущения следует учитывать при определении орбиты из наблюдений и при прогнозировании движения опасного тела? Совершенно очевидно, что следует учитывать гравитационные возмущения тела от больших планет Солнечной системы и Луны. Для получения координат и скоростей планет, Земли и Луны в настоящее время повсеместно используются численные эфемериды этих тел, публикуемые Лабораторией реактивного движения (ЛРД) НАСА (DE405/LE405 [Standish, 2000] и более поздние промежуточные версии). Численные эфемериды того же уровня точности EPM2004, EPM2006, EPM2008 созданы и постоянно совершенствуются в Институте прикладной астрономии РАН [Питьева, 2005; 2007]. В основе всех этих численных теорий лежат релятивистские (составленные на основе общей теории относительности) уравнения движения больших планет, Солнца и Луны относительно барицентра Солнечной системы [Расширенное объяснение к Астрономическому ежегоднику, гл. 6, 2005]. Многочисленные параметры теории движения (их около 200) определяются из высокоточных радарных наблюдений внутренних планет, радионаблюдений КА, лазерных наблюдений Луны и оптических наблюдений планет и спутников, полученных за последние 100 лет. В уравнения движения тел включены возмущения, вызываемые сжатием Солнца, и возмущения от 300 наиболее массивных астероидов. Найденные значения параметров теории прошли многочисленные перекрестные проверки и сравнения. Точность теории DE405 оценивается величиной порядка 1 км для внутренних планет. Надо отметить, что для разных планет и Луны численные эфемериды обеспечивают различную точность. Некоторое представление о точности координат планет, вычисляемых по этим эфемеридам, можно составить, сравнивая между собой максимальные различия значений гелиоцентрических расстояний планет, вычисляемых по теориям, созданным в разных организациях. Так, например, на интервале с 1960 г. по 2020 г. модуль разности этих величин, вычисляемых по DE405 и EPM2006, для Венеры не превышает 180 м, для Земли — 26 м, для Марса — 120 м, для Юпитера — 16 км [Питьева, 2007]. Более поздние версии эфемерид обнаруживают еще меньшие различия.
Вопрос о том, насколько сказывается неточность используемых координат возмущающих планет, Земли и Луны на прогнозе столкновений опасных космических тел с Землей, может быть исследован на примере все того же опасного астероида Апофис. С этой целью следует провести уточнение параметров орбиты астероида и последующее прогнозирование его движения дважды: один раз с использованием, например, численной эфемериды DE405, а второй раз с использованием более современной эфемериды DE414. Такие вычисления были выполнены в ЛРД [Giorgini et al., 2008] и в ИПА РАН на основе несколько различающихся моделей движения. Из результатов этих вычислений следует, что различие между DE405 и DE414 оказывает заметное влияние на координаты Апофиса в апреле 2029 г. Разница расстояний до центра Земли, рассчитанная с применением этих двух теорий, на входе в сферу действия Земли составляет всего –0,597 км. Та же величина, рассчитанная в момент наибольшего сближения тел с Землей, составила уже +3,2 км. Это различие продолжает стремительно возрастать при переходе к сближению в апреле 2036 г., когда оно уже составляет 136 000 км. Причина заключается в том, что в ходе тесного сближения в апреле 2029 г. Земля оказывает очень разное воздействие на тела, движущиеся по орбитам с минимальными геоцентрическими расстояниями, различающимися всего на несколько километров. Пучок орбит с почти одинаковыми геоцентрическими расстояниями после тесного сближения расходится своеобразным веером, вследствие чего точность последующего предсказания резко снижается. Это весьма неприятная особенность тесных сближений, поскольку близкие прохождения астероидов около Земли чреваты повторными сближениями с ней спустя всего лишь несколько лет (см. раздел 7.7.1), и надо уметь заранее предвычислять эти сближения. Потеря точности при тесных сближениях предъявляет весьма суровые требования к точности исходной орбиты.
Выше было отмечено, что теории DE и EPM в настоящее время строятся с учетом возмущений от 300 малых планет. Однако учет этих возмущений оказывает незначительное влияние на движение Земли, Луны и других планет. Поэтому при исправлении орбит подавляющего числа потенциально опасных астероидов и при прогнозе их движения целесообразно включать в число возмущающих планет только большие планеты (возмущения от Земли и Луны учитываются раздельно) и три наиболее массивных астероида: Цереру (которая в настоящее время рассматривается как карликовая планета), Палладу и Весту. Точность, с которой вычисляются координаты трех последних тел, не является при этом критически важной. Величины возмущающих масс, обязательные для использования для всех возмущающих тел, указываются в описаниях соответствующих теорий.
Влияние несферичности гравитационного поля Земли на движение Апофиса до и после сближения 2029 г. исследовалось численным образом в работе [Giorgini et al., 2008]. При отсутствии тесных сближений с Землей влияние ее несферичности на движение астероида является минимальным. Однако прогноз сближения в 2036 г. при использовании сферической модели Земли дает ошибку порядка 19 000 км по сравнению с прогнозом, учитывающим ее сжатие. Прочие зональные гармоники вносят возмущения, по крайней мере, на три порядка меньшие. Долготная вариация гравитационного поля Земли приводит к уклонению порядка 100 км в 2036 г.
Влияние возмущений, вызываемых тепловым излучением астероида (эффект Ярковского), рассмотрено в разделе 7.7.6.
Обратимся теперь к проблеме численного интегрирования уравнений движения малых тел. Существует большое число методов численного интегрирования дифференциальных уравнений движения, способных обеспечить достаточно высокую точность прогноза движения опасного астероида на длительных интервалах времени (см., напр., [Бордовицына, 1984; Бордовицына, Авдюшев, 2007]). В различных организациях для этой цели используются методы Булирша — Штера, Эверхарта, Рунге — Кутты — Фельберга, Адамса и другие. В ЛРД НАСА на протяжении многих лет широко используется метод с переменным шагом Адамса — Крога для интегрирования дифференциальных уравнений второго порядка [Krogh, 1968; 1974]. Метод применяется при построении численных эфемерид больших планет и Луны при проведении многочисленных космических экспедиций к телам Солнечной системы. Он многократно доказал свою точность и эффективность. Метод реализован с двукратной точностью (ошибка порядка 10-14 на каждом шаге) и даже четырехкратной точностью (ошибка порядка 10-19 на одном шаге). Правда, четырехкратная точность требует примерно в 30 раз больше компьютерного времени.
В ИПА РАН при построении численных эфемерид больших планет и Луны, как и при расчетах эфемерид малых тел, широко используется метод численного интегрирования Эверхарта [Everhart, 1974a; 1974b], интегратор RADAU с двойной точностью. Метод позволяет, в частности, интегрировать уравнения вида (зависимость правых частей от появляется при учете релятивистских членов и негравитационных эффектов). Метод дает возможность вести вычисления с различной точностью в зависимости от учитываемого числа членов в разложениях и числа последовательных приближений на одном шаге.
Интересно получить ответ на вопрос, какая часть ошибки прогноза на 2029 г. и на 2036 г. может быть связана с ошибками численного интегрирования. Представление об этом дают «ошибки замыкания», вычисляемые при интегрировании уравнений движения Апофиса на различных интервалах времени при задании различных порядков интегрирования. Под ошибкой замыкания понимается модуль разности между исходным значением гелиоцентрического радиуса-вектора астероида и его окончательным значением в тот же самый момент, если интегрирование выполняется от исходной эпохи до конечной, а затем назад к исходной. Исследование было выполнено на трех интервалах: 1) от 2007 г. до 12 апреля 2029 г. (до входа Апофиса в сферу действия Земли); 2) от 2007 г. до момента наибольшего сближения Апофиса с Землей в 2029 г.; 3) от 2007 г. до 2036 г. Оказалось, что выбор порядка интегрирования слабо влияет на достигаемую точность. На первом интервале ошибка замыкания не превышает или равна 10-11 а.е., что вполне достаточно для сравнения вычисляемых положений с наблюдаемыми не только для оптических, но и для радиолокационных наблюдений. Погрешность вычисления геоцентрического расстояния астероида, зависящая от численного интегрирования, на границе сферы действия (до сближения в 2029 г.) составляет 2 × 10-5 км. Точность вычислений заметно снижается при тесном сближении (второй интервал). Ошибка замыкания возрастает до 10-8 а.е., но минимальное расстояние от Земли вычисляется с точностью до 30 см (ошибка, зависящая только от численного интегрирования). Наконец, к концу третьего интервала после сближения в 2029 г. координаты астероида вычисляются с ошибками порядка нескольких километров, а ошибка замыкания на этом интервале возрастает до 10-7 а.е. Таким образом, численное интегрирование не вносит существенных ошибок в вычисляемые положения на рассматриваемых интервалах.
7.4. Траектория сближения тела с Землей и другими массивными телами. Гравитационный маневр. Радиус захвата. Плоскость цели
При оценке вероятности столкновения естественных космических тел друг с другом или искусственных космических аппаратов с естественными телами важнейшую роль играет понятие плоскости цели. Плоскость цели — это плоскость, проходящая через центр планеты-мишени перпендикулярно к вектору невозмущенной скорости тела-снаряда относительно планеты-мишени. Когда астероид имеет тесное сближение с большой планетой, его гелиоцентрическая орбита начинает постепенно меняться под действием тяготения планеты. Внутри сферы действия планеты траектория астероида относительно планеты очень близка к гиперболе (рис. 7.1) (напомним, что сферой действия планеты называется область пространства, в которой отношение возмущающего ускорения, сообщаемого телу планетой, к ускорению, сообщаемому телу Солнцем, превосходит отношение возмущающего ускорения, сообщаемого телу Солнцем, к ускорению, сообщаемому телу планетой; приближенное значение радиуса сферы действия Земли равно 0,0062 а.е., или 930 000 км).
Скорость астероида относительно Земли на входе в сферу действия на разности гелиоцентрических скоростей астероида и Земли. Это так называемая скорость тела относительно Земли на бесконечности (невозмущенная скорость тела относительно Земли). По направлению она близка к асимптоте гиперболы, описываемой телом в сфере действия планеты (рис. 7.1).
Рис. 7.1. Траектория движения астероида относительно Земли в пределах ее сферы действия
Обогнув Землю (как говорят, совершив гравитационный маневр), на выходе из сферы действия астероид имеет ту же самую по величине относительную скорость , но ее направление изменяется на угол γ. Гелиоцентрическая скорость тела на выходе из сферы действия в результате поворота вектора скорости также меняется.
Из определения плоскости цели следует, что на рис. 7.1 штриховая прямая, проведенная перпендикулярно асимптоте гиперболы через центр Земли, есть след от пересечения плоскости цели с плоскостью орбиты тела относительно Земли. Отрезок этой прямой от центра Земли до асимптоты обозначен как b. Его называют прицельным расстоянием.
Как видно из рисунка, прицельное расстояние по величине превышает минимальное расстояние от гиперболы до центра Земли q. Эти две величины связаны соотношением
где v∞ есть параболическая скорость относительно Земли:
Здесь G — гравитационная постоянная, M⊕ — масса Земли, r⊕ — ее экваториальный радиус. Если в формулу (7.9) подставить q, равное r⊕, то b будет равно прицельному расстоянию, при котором траектория астероида коснется поверхности Земли. Соответствующее значение прицельного расстояния называется радиусом захвата. При меньших значениях прицельного расстояния астероид обязательно столкнется с Землей. В зависимости от соотношения и v∞ радиус захвата может существенно превышать геометрический радиус Земли. При решении вопроса о реальности столкновения следует в некоторых случаях использовать не радиус Земли, а ее радиус захвата.
Рассмотрение процесса сближения космических тел с Землей облегчается при использовании специально выбранной системы координат. Столкновения могут иметь место только в малой окрестности минимального расстояния между орбитами. В этой окрестности орбиты Земли и тела могут рассматриваться как отрезки двух прямых, скрещивающихся в пространстве (в частном случае — пересекающихся). Кратчайшим расстоянием между ними является отрезок прямой, перпендикулярный к обеим скрещивающимся прямым.
При выборе системы координат ее начало помещают в центр Земли. Плоскость цели проводят через центр Земли перпендикулярно к вектору геоцентрической скорости астероида . Отметим особо два момента. Первый момент: отрезок кратчайшего расстояния между двумя орбитами EA лежит в плоскости цели (рис. 7.2). Второй момент: в малой окрестности кратчайшего расстояния между орбитами, где орбиты могут рассматриваться как отрезки прямых, движение обоих тел происходит в параллельных плоскостях.
На рис. 7.2 EE1 — орбита Земли, AA1 — орбита астероида, EA — отрезок кратчайшего расстояния между двумя орбитами; V — гелиоцентрическая скорость Земли в тот момент, когда она проходит через точку E, v — гелиоцентрическая скорость виртуального астероида, который проходит через точку A в тот момент, когда Земля оказывается в точке E. Астероид, соответствующий номинальному решению, проходит через точку A, вообще говоря, раньше или позже того момента, когда Земля находится в точке E, но эти моменты времени не сильно отличаются друг от друга, иначе столкновение не происходит.
Рис. 7.2. Система координат, связанная с плоскостью цели. Показаны положения плоскости цели в два разных достаточно близких момента времени
На рисунке показана также скорость виртуального астероида относительно Земли . Ось η системы координат направлена параллельно , ось ξ направлена по векторному произведению скорости Земли V на орт оси η. Из этого следует, что ось ξ направлена вдоль отрезка EA кратчайшего расстояния между орбитами. Ось ζ направлена так, чтобы она вместе с осями ξ и η образовывала правую систему координат. Угол θ на рисунке есть угол в плоскости η — ζ между направлением скорости Земли и осью η.
Формулы перехода от системы координат x, y, z к системе ξ, η, ζ имеют вид
где x0, y0, z0 — координаты центра Земли.
Направляющие косинусы осей ξ, η, ζ относительно x, y, z находятся из соотношений
В точке A орбита выбранного нами виртуального астероида пересекает плоскость цели. Будем вести отсчет времени от этого момента пересечения. Пусть выбранный астероид движется по орбите, соответствующей номинальному решению, но значение средней аномалии для него отлично от значения средней аномалии в номинальном решении. Расстояние EA, или координата ξ точки A, равны наименьшему возможному расстоянию орбиты от центра Земли. В момент пересечения плоскости цели координата ζ равна нулю. Но с течением времени координата ζ точки пересечения данного астероида с плоскостью цели не остается постоянной, поскольку начало системы координат движется вместе с Землей. По истечении промежутка времени Δt координата ζ точки, в которой некоторое время тому назад произошло пересечение астероида с плоскостью цели, будет равна
ζ = |V | sin θ Δt, (7.11)
где θ — угол между направлением гелиоцентрической скорости Земли и осью η.
В отличие от координаты ζ координата ξ не меняется с течением времени, поскольку движение происходит в параллельных плоскостях. Если пересечение с плоскостью цели рассматривать на границе сферы действия планеты, то координата ξ равна прицельному расстоянию b (рис. 7.1). Если пересечение рассматривается внутри сферы действия планеты вблизи перигея гиперболы, то ξ равно минимальному расстоянию q от гиперболы до центра Земли.
7.5. Эллипс рассеяния в плоскости цели. Оценка вероятности столкновения
Только один виртуальный астероид пересекает плоскость цели в момент, когда Земля находится у одного конца кратчайшего отрезка между орбитами. Другие виртуальные астероиды, движущиеся вдоль номинальной траектории, пересекают плоскость цели раньше или позже, чем это нужно для достижения минимального расстояния между орбитами, и соответствующие точки пересечения имеют различные значения координаты ζ. Очевидно, что есть расстояние, на котором виртуальный астероид пересекает плоскость цели от центра Земли. В то же время это есть минимальное расстояние от Земли, на котором он проходит мимо нее в данном сближении. Это следует из того, что его геоцентрическая скорость нормальна к геоцентрическому радиусу.
Таким образом, цепочка виртуальных астероидов, вытянувшихся вдоль номинальной орбиты, проектируется на плоскость цели в прямую, параллельную оси ζ, причем виртуальный астероид, соответствующий центру доверительного эллипсоида в начальную эпоху t0, пересекает плоскость цели в точке, расположенной, вообще говоря, выше или ниже оси ξ. Область вокруг этой точки на плоскости ξ — ζ является отображением области возможных начальных условий движения на плоскость цели. Поскольку мы с самого начала предположили линейный характер задачи, можно утверждать, что область начальных значений, ограниченная в эпоху t0 доверительным эллипсоидом, отобразится на плоскость ξ — ζ в часть плоскости, ограниченную эллипсом с центром в точке, соответствующей центру доверительного эллипсоида. Задача сводится к тому, чтобы найти координаты центра эллипса на плоскости ξ — ζ и его полуоси и оценить расположение эллипса рассеяния относительно образа Земли на этой плоскости.
В линейном приближении эта задача решается достаточно просто. В общем виде ход решения задачи можно описать следующим образом.
Координаты точки ξ, ζ на плоскости цели (см. формулу (7.10)) являются функциями F1 и F2 параметров орбиты (элементов или координат и скоростей в начальную эпоху), что в векторном виде можно записать как
L = F(E),
где L — двумерный вектор с компонентами ξ, ζ, а E — вектор параметров орбиты.
В рамках линейного приближения матрица ковариации D вектора L связана с матрицей ковариации вектора E известным соотношением [Эльясберг, 1976]:
D = σ 2 ( d F d E)Q -1 ( d F d E) T , (7.12)
где dF/dE — частные производные F по параметрам E. Величины σ и Q-1 известны, поскольку они являются соответственно средней квадратичной погрешностью наблюдений, использованных при определении орбиты тела, и обратной матрицей нормальной системы уравнений (см. раздел 7.1).
Компоненты вектора L и частные производные в момент t (изохронные производные) находятся численным интегрированием уравнений движения в прямоугольных координатах с последующим преобразованием их в координаты ξ, η, ζ и численным интегрированием уравнений, определяющих значения производных (так называемых уравнений в вариациях) при заданных начальных условиях движения. Таким образом, на момент сближения астероида, соответствующего номинальному решению, с Землей (или со сферой ее действия) оказываются известными координаты центра эллипса в плоскости цели и его полуоси, определяемые как
где Dii — диагональные элементы матрицы ковариации D, a1 = aξ — длина малой полуоси эллипса рассеяния, a2 = aζ — длина большой полуоси. Заметим, что формула (7.13) определяет полуоси эллипса, соответствующие области неопределенности начальных условий внутри эллипсоида равных плотностей вероятности. Чтобы получить полуоси доверительного эллипса на плоскости цели, надо aξ и aζ умножить на 3.
Возможны следующие три случая взаимного расположения Земли и эллипса на плоскости цели: а) эллипс расположен на некотором расстоянии от окружности с радиусом, равным радиусу Земли (радиусу захвата Земли, если вычисления доверительного эллипса производятся на границе сферы действия Земли) (рис. 7.3 а), что практически исключает возможность столкновения астероида с Землей;
б) кружок с радиусом, равным радиусу Земли (или радиусу захвата), находится внутри эллипса (рис. 7.3 б). Вероятность столкновения может быть рассчитана исходя из отношения площади кружка к площади, ограниченной эллипсом. Для повышения точности прогноза можно учесть неодинаковую вероятность попадания виртуальных астероидов в различные точки области, ограниченной эллипсом;
Рис. 7.3. Возможные взаимные расположения эллипсов рассеяния и Земли в плоскости цели
в) площадь, ограниченная эллипсом, частично покрывает Землю (рис. 7.3 в). Этот случай практически не отличается от предыдущего. Вероятность столкновения рассчитывается с учетом отношения перекрывающейся области ко всей площади, ограниченной эллипсом.
Более подробно расчет вероятности столкновения здесь не рассматривается, так как во всех случаях, когда возникает реальная угроза столкновения, следует предпринять дополнительные исследования, учитывающие возможный нелинейный характер задачи.
Нелинейный характер задача может иметь по многим причинам. Доверительный эллипсоид уже в эпоху t0 может недостаточно хорошо описывать область возможных начальных условий, поскольку само распределение ошибок наблюдений может не подчиняться закону Гаусса. Чем дальше от эпохи t0, тем больше нарастает нелинейность, и применение формулы (7.8) становится незаконным. Проекция доверительного эллипсоида на плоскость цели в момент t сближения с Землей, отдаленный от t0 на десятилетия, вытягивается в очень узкую область, которая к тому же искривляется в соответствии с кривизной земной орбиты. По всем этим причинам линейный анализ задачи становится неадекватным и требуется применение более тонких методов анализа. К настоящему времени предложено два таких метода: метод Монте-Карло и метод линии вариации.
Метод Монте-Карло, или метод статистических испытаний, в применении к данной задаче означает прямое использование вероятностной интерпретации метода наименьших квадратов. Поскольку процесс уточнения орбиты по МНК доставляет, как принято говорить, наиболее вероятное решение, окруженное областью других возможных решений, то можно выбрать в этой области случайным образом большое число виртуальных астероидов и следить со всей возможной точностью за их движением в течение некоторого времени, пока они не столкнутся с Землей или не пролетят мимо нее. Тогда отношение числа столкнувшихся виртуальных астероидов к их общему количеству можно рассматривать как вероятность столкновения с Землей астероида, орбита которого доподлинно неизвестна. Этот метод замечателен своей простотой, универсальной применимостью и правильным учетом нелинейности задачи. При его практическом использовании важно учитывать корреляционные зависимости между разыгрываемыми значениями параметров орбиты, но это реализуется достаточно просто [Железнов, 2009]. К сожалению, метод является чрезвычайно трудоемким. Действительно, чем менее вероятное событие требуется оценить, тем большее количество начальных условий движения следует испытать. Пусть, например, при испытании 106 случайно выбранных начальных условий в пяти случаях было зафиксировано столкновение с Землей. Тогда можно утверждать, что вероятность столкновения близка к 0,000005. Но если проведена только тысяча испытаний, которые не дали ни одного попадания, тогда можно лишь сказать, что вероятность столкновения, по-видимому, меньше 0,001. Поскольку на практике приходится искать опасные сближения с Землей на интервалах в несколько десятков лет и вероятность столкновения при этом имеет, как правило, порядок 10-4 и менее, то требуется несколько дней работы компьютера для получения надежного результата в отношении только одного астероида [Milani et al., 2000].
Метод Монте-Карло основывается на выборе случайных точек во всем шестимерном пространстве возможных начальных условий и на их последующем испытании. Имеется также возможность выбора точек в каком-нибудь подпространстве, относительно которого можно предполагать, что берущие в нем начало решения достаточно хорошо отражают поведение решений во всей доверительной области. В качестве такого подпространства можно, например, использовать линию вариации, вдоль которой номинальное решение определяется с наибольшей погрешностью. В доверительном эллипсоиде линия вариации совпадает с направлением наиболее вытянутой оси, как правило, большой полуоси его орбиты. В методе линии вариации виртуальные астероиды берутся со значениями пяти элементов, соответствующими номинальному решению, в то время как шестой элемент (среднее движение или большая полуось) варьируется с постоянным шагом в пределах ±3σ (или в иных пределах). Как и в методе Монте-Карло, движение виртуальных астероидов прослеживается на всем исследуемом интервале, в особенности при их сближениях с Землей. Поскольку при этом точки пересечения виртуальных астероидов с плоскостью цели представляют наборы, зависящие только от одного параметра, то достаточно просто (путем интерполяции или методом Ньютона нахождения корней функции) определяются значения среднего движения (большой полуоси), при которых реализуется максимальное сближение виртуального астероида с Землей.
Хотя этот метод является эффективным средством анализа сближений, нельзя быть уверенным, что при этом будут найдены все возможные столкновения, например те, которые соответствуют точкам доверительного эллипсоида, расположенным далеко от линии вариации. Соответствующие им точки на плоскости цели, если имеет место сильно выраженная нелинейность задачи, могут оказаться на значительном удалении от точек, отвечающих линии вариации, и часть из них может при этом вести к столкновениям. Метод Монте-Карло должен, в принципе, обнаруживать подобные случаи. Поэтому оба метода должны дополнять друг друга и использоваться для взаимного контроля.
7.6. Потоки виртуальных астероидов, следующие различными динамическими путями
Обратимся теперь к рассмотрению особенностей, которые связаны с нелинейными эффектами. Если все возможные сближения виртуальных астероидов с Землей на исследуемом интервале упорядочить по времени, то можно видеть, что они группируются около нескольких эпох, когда Земля оказывается вблизи узла номинальной орбиты астероида на эклиптике. По аналогии с метеороидами, встречающимися с Землей, такие наборы виртуальных астероидов можно назвать потоками. Очень часто поток можно подразделить на отдельные струи, охватывающие подмножества виртуальных астероидов на динамически различных путях. Например, один поток может включать несколько струй, в которых виртуальные астероиды совершили различное число оборотов вокруг Солнца за время, истекшее с момента t0. Аналогичное подразделение потока на струи может возникнуть в результате его тесного сближения с Землей, когда часть виртуальных астероидов под влиянием возмущений со стороны последней оказывается на орбитах, имеющих те или иные соизмеримости с Землей. В итоге, спустя определенное целое число лет, Земля может встретиться вблизи узла номинальной орбиты астероида уже с несколькими различными струями виртуальных астероидов. На плоскости цели отдельные струи обнаруживаются в виде чрезвычайно вытянутых цепочек виртуальных астероидов, располагающихся внутри доверительной области отдельной струи, ширина которой во много раз меньше ее длины.
Поток виртуальных астероидов в январе 2046 г., соответствующих орбите астероида 1998 OX4, подразделяется на три струи. Одна из струй пересекает Землю, ввиду чего возможно столкновение. Некоторые струи могут обнаруживать поведение, отличное от только что описанного. Например, при изменении среднего движения вдоль линии вариации точки пересечения соответствующих виртуальных астероидов с плоскостью цели сначала приближаются к линии минимального расстояния между Землей и астероидом, но затем движение прекращается и сменяется на обратное. Такое поведение может быть названо прерванным возвращением. Оно возникает в результате предшествующих тесных сближений с Землей. Примером прерванного возвращения может служить струя потока виртуальных астероидов в январе 2046 г., соответствующих орбите астероида 1998 OX4 [Milani et al., 2002].
Остановимся на некоторых особенностях, связанных с подсчетом вероятности столкновения в случае нелинейности задачи. Как мы видели, для таких случаев характерно вытягивание области пересечения виртуальных астероидов с плоскостью цели в очень узкую полосу, длина которой в тысячи раз превосходит ширину. Точка, соответствующая номинальному решению, может отстоять на очень большое расстояние от Земли. При этом из-за нелинейности точность определения минимального расстояния полосы от Земли оказывается невысокой, что недопустимо, если само расстояние составляет, согласно вычислениям, всего лишь несколько земных радиусов. В таких случаях рекомендуется использовать итеративную процедуру, которая сводится к поиску поправок к номинальному решению, приводящих к более тесному сближению соответствующего виртуального астероида с Землей [Milani et al., 2002]. Процесс повторяют до тех пор, пока минимальное расстояние не перестанет изменяться.
Если оказалось, что минимально возможное расстояние меньше радиуса захвата Земли, то столкновение в принципе возможно и надо оценить его вероятность. При условии, что ширина полосы много меньше диаметра Земли, легко написать формулу для оценки вероятности столкновения. В самом деле, полоса, где располагаются виртуальные астероиды, вырезает в области захвата Земли хорду l максимальной длины 2R0 (R0 — радиус захвата Земли). Пусть длина всей полосы, соответствующая изменению среднего движения в пределах ±3σn составляет Λ км (σn — средняя ошибка среднего движения). Можно написать, что Λ = 6aζ, где aζ — длина большой полуоси эллипса рассеяния (см. формулу (7.13)).
Плотность вероятности распределения виртуальных астероидов вдоль этой полосы, которую мы будем рассматривать как одномерную, зависит от двух параметров: центра распределения и величины aζ. Центр распределения в данном случае совпадает с серединой полосы (точкой пересечения с плоскостью цели того виртуального астероида, который имеет среднее движение, найденное в номинальном решении). Величина aζ также известна. Плотность распределения виртуальных астероидов вдоль полосы описывается кривой Гаусса
где p(z) — плотность вероятности при данном значении z, отсчитанном от z0. Величина p(z)dz есть вероятность попадания виртуального астероида в интервал z — (z + dz). Чтобы воспользоваться этой формулой, надо оценить величину p(z) при значении z, соответствующем пересечению полосы с Землей, и величину dz. Отношение l/aζ, где l — длина хорды, можно рассматривать как малое по величине приращение dz переменной z в пределах области захвата. Значение переменной z, соответствующее пересечению полосы с Землей, — это просто расстояние от центра полосы до центра Земли, выраженное в единицах aζ.
В случае, если ширина полосы существенно больше диаметра Земли, для оценки вероятности столкновения приходится уже решать двумерную задачу с учетом неравномерной плотности вероятности как вдоль полосы, соответствующей изменению среднего движения вдоль линии вариации, так и поперек нее. Формула для вычисления вероятности столкновения в этом случае приобретает вид
где p1 и p2 — плотности вероятности вдоль полосы и поперек нее в точке, соответствующей центру Земли, вычисляемые аналогично (7.14), aζ и aξ — соответственно длина большой и малой полуоси эллипса рассеяния, а R0 — радиус захвата Земли.
Прогнозы возможных столкновений и ссылки на сайты соответствующей тематики приведены в приложениях 1, 3, 4.
7.7. Опасный астероид (99942) Апофис
7.7.1. История обнаружения и исследования астероида Апофис. Несколько следующих подразделов посвящены результатам исследования движения потенциально опасного астероида Апофис с применением описанных выше методов.
Данный астероид, получивший предварительное обозначение 2004 MN4, был открыт в обсерватории Китт Пик (штат Аризона, США) 19 июня 2004 г. в ходе регулярных наблюдений, проводившихся там по программе поиска тел, сближающихся с Землей. Было получено шесть наблюдений астероида в течение двух ночей, после чего он фактически был потерян. 18 декабря того же года этот астероид был случайно переоткрыт Г. Гарраддом в Австралии, выполнявшим наблюдения по программе обзора астероидов в обсерватории Сайдинг Спринг. На этот раз наблюдения были продолжены во многих обсерваториях, что дало Центру малых планет возможность уже 20 декабря опубликовать Электронный циркуляр c орбитой астероида, найденной на основе 31 наблюдения, выполненного с июня по декабрь. 2004 MN4 оказался астероидом, сближающимся с Землей, движущимся по орбите, внутренней по отношению к земной и мало наклоненной к ней. Минимальное расстояние его орбиты от орбиты Земли составило всего 0,002 а.е. Абсолютная звездная величина астероида была оценена как 19,3m, что при характерных для данного типа астероидов альбедо соответствовало диаметру около 400 м. Малость межорбитального расстояния и размеры астероида дали основание классифицировать его как потенциально опасный для Земли.
В течение ближайших трех дней после повторного открытия астероида сразу две автоматические системы SENTRY в США и CLOMON 2 в Италии нашли, что астероид будет иметь очень тесное сближение с Землей 13 апреля 2029 г. Первоначальные оценки вероятности его столкновения с Землей в 2029 г. составляли около одной трехсотой, причем по мере увеличения числа используемых для уточнения орбиты наблюдений вероятность столкновения только возрастала. К 26 декабря по 169 наблюдениям вероятность оценивалась уже как 2,2 %, а к 27 декабря по 176 наблюдениям — как 3,3 %. Напряженность ситуации нарастала. Столкновение Земли с астероидом таких размеров угрожало региональной катастрофой.
К счастью, к вечеру 27 декабря в обсерватории Китт Пик были обнаружены шесть слабых изображений астероида 2004 MN4, полученных по иной программе наблюдений в марте 2004 г. и ранее остававшихся незамеченными. Расширение интервала времени, охваченного наблюдениями, позволило значительно увеличить точность предсказаний и по существу исключить возможность столкновения Земли с данным астероидом в 2029 г. Тем не менее, астероид продолжал интенсивно наблюдаться. В конце января 2005 г. в радиообсерватории Аресибо (Пуэрто-Рико) была проведена серия из пяти радиолокационных наблюдений астероида, позволивших существенно уточнить элементы его орбиты и с большей достоверностью прогнозировать его движение в будущем. Еще два радарных наблюдения астероида были выполнены 7 августа 2005 г. и 6 мая 2006 г. Всего на интервале с 15 марта 2004 г. по 16 августа 2006 г., по данным Центра малых планет, было выполнено 996 наблюдений, в том числе 989 наблюдений прямого восхождения и склонения и семь радарных наблюдений: два наблюдения запаздывания отраженного сигнала и пять наблюдений доплеровского сдвига частоты (http://ssd.jpl.nasa.gov/?radar).
Как оказалось, астероид пройдет мимо Земли 13 апреля 2029 г. на расстоянии даже меньшем, чем это ранее предполагалось: от центра Земли оно составит немного более 38 000 км — меньше высоты полета геостационарных спутников. На ночной половине Земли астероид можно будет наблюдать как звезду 3,3m. Для землян столь близкое прохождение небесного тела само по себе мало чем опасно. Однако опасными могут оказаться его последствия. Для астероида сближение с Землей обернется большим или меньшим изменением элементов орбиты под влиянием возмущения со стороны гравитационного поля планеты. Величина возмущений существенно зависит от минимального расстояния до астероида. К сожалению, сегодня это расстояние может быть рассчитано недостаточно точно. Отсюда возникает некоторая неопределенность прогноза дальнейшего движения астероида. В любом случае значительно изменится большая полуось его орбиты, или, другими словами, его среднее движение. Орбита астероида из внутренней по отношению к Земле (орбиты типа Атона, выходящей за пределы орбиты Земли только в окрестности афелия) превратится в орбиту типа Аполлона, с большой полуосью, превышающей 1 а.е., и подходящую к Солнцу ближе орбиты Земли в окрестности своего перигелия. Орбитальный период, который в настоящее время равен 0,89 года, превысит величину одного года. При этом может оказаться, что продолжительность нового периода, выраженная в годах, составит величину, представляющую собой отношение двух целых чисел, например 7: 6. Это будет означать, что по прошествии семи лет тело совершит относительно Солнца шесть оборотов, в то время как Земля совершит в точности семь оборотов. По истечении этого периода времени ситуация встречи тел на орбитах должна повториться. При небольшом отличии от точной соизмеримости периодов возможно даже ухудшение ситуации, имевшей место в предыдущем сближении. Например, если в предыдущем сближении тело прошло через точку минимального расстояния между орбитами немного ранее, чем Земля, ввиду чего столкновения тел не произошло, то к моменту очередного сближения тело может несколько задержаться и прийти в роковую точку одновременно с Землей. Как показывают расчеты, для 2004 MN4 в результате сближения в 2029 г. возможен целый спектр неблагоприятных преобразований орбиты, которые ведут либо к непосредственным столкновениям через несколько лет с Землей, либо к новым очень тесным сближениям с нею и новым опасным трансформациям орбиты в будущем.
В ближайшие десятилетия или даже столетия астероид будет являться источником постоянной угрозы для Земли, пока вековые возмущения перигелия и узлов орбиты не уведут его восходящий узел на эклиптике достаточно далеко от орбиты Земли. Правда, эти же возмущения через несколько тысяч лет приведут к опасной близости к орбите Земли его нисходящий узел [Заботин, Медведев, 2008]. На интервалах в десятки тысяч лет движения перицентра и узлов будут периодически приводить орбиту астероида к пересечениям с орбитой Земли. Орбита астероида является, таким образом, динамически неустойчивой. В конечном итоге астероид либо выпадет на Землю, либо в результате тесных сближений с нею большая полуось и эксцентриситет его орбиты будут радикально преобразованы, и астероид выпадет на Солнце (при очень большом эксцентриситете) или получит возможность сближаться с другими большими планетами и эволюционировать под влиянием тесных сближений с ними.
Такова в общих чертах картина движения данного астероида и эволюции его орбиты. В июне 2005 г. астероид 2004 MN4 получил постоянный номер 99942 и вскоре после этого ему было присвоено имя Apophis (Разрушитель) (греческое имя египетского божества Aпеп, которое, согласно легенде, таится во мгле и пытается разрушить Солнце во время его ночных странствий по подземному миру).
7.7.2. Точность знания орбиты Апофиса на современном этапе. Для оценки угрозы Земле со стороны (99942) Апофиса первостепенное значение имеет знание его точной орбиты. Орбита астероида имеет необычайно малое межорбитальное расстояние с орбитой Земли, равное 0,00027 а.е. Это определяет возможность очень тесных сближений астероида с Землей, которые происходят в окрестности 13 апреля, если Земля пересекает линию узлов орбиты Апофиса на эклиптике почти одновременно с астероидом. Ближайшее по времени тесное сближение Апофиса с Землей произойдет 13 апреля 2029 г., когда минимальное расстояние между их центрами будет около 38 000 км (рис. 7.4). Размеры Апофиса (∼270 ± 60 м) и особенности его орбиты таковы, что до 2012 г. трудно ожидать его новых наземных наблюдений (см. ниже).
Рис. 7.4. Положения Земли и Апофиса на их орбитах 1 апреля 2029 г.
Радарные наблюдения более информативны и позволяют даже при небольшом их числе существенно повысить точность определения орбиты. Известны, по крайней мере, четыре определения орбиты Апофиса с использованием имеющихся радарных и большей части оптических наблюдений. Все четыре системы параметров орбиты найдены путем взвешенного уравнивания систем условных уравнений и их решения по методу наименьших квадратов. Представляет значительный интерес сравнение найденных параметров орбит и полученных оценок их точности. Как уже было сказано, из-за небольших различий в значениях параметров в начальную эпоху проистекают существенные различия в прогнозе движения Апофиса после 2029 г. и различия в оценке вероятности его столкновения с Землей в 2036 г. и в последующие годы. Ниже дается оценка того, насколько близки найденные решения друг к другу и сколь существенно может отличаться прогноз движения на время после 2029 г. из-за различия моделей движения и неучтенных факторов, влияющих на движение астероида.
Были проведены сравнения систем элементов орбиты и оценки их точности, полученные в Институте прикладной астрономии (ИПА) РАН, Лаборатории реактивного движения (ЛРД НАСА, США), в Пизанском университете (Италия) (NEODyS) и группой сотрудников ЛРД, радиообсерватории Аресибо и Калифорнийского технологического института [Giorgini et al., 2008]. Первые три системы найдены на эпоху 2007, апрель 10,0, четвертая система — на эпоху 2006, сентябрь 1,0. При получении каждой системы было использовано разное число оптических наблюдений. Различалась и процедура назначения весов условным уравнениям. Соответственно различаются и оценки средней ошибки. Тем не менее, первые три системы весьма близки друг к другу, и различие между ними находится в пределах найденных ошибок элементов.
В табл. 7.1 представлены значения координат и компонент скорости Апофиса и их среднеквадратичные погрешности (решение ИПА РАН). Для сравнения рядом указываются данные, соответствующие решению, полученному только по оптическим наблюдением. Очевидно, что включение в решение радарных наблюдений увеличивает его точность в несколько раз.
Таблица 7.1. Координаты и компоненты скорости Апофиса в эпоху JD 2454200,5 (10.04.2007) и их средние ошибки, полученные при использовании оптических и радарных и только оптических наблюдений
7.7.3. Сближение с Землей в 2029 г. Интересно проследить, какие следствия вытекают из различий между полученными решениями. Эти различия в первую очередь касаются минимального расстояния между Землей и Апофисом, достигаемого 13 апреля 2029 г. Величины этих расстояний в момент времени JD 2462240,407115, близкий к эпохе максимального сближения этих тел в 2029 г., представлены в табл. 7.2.
Следует оговориться, что только в случае решения ИПА эти величины подсчитаны в полном соответствии с моделью, использованной при уточнении орбиты. В решении [Giorgini et al., 2008]) использована так называемая стандартная динамическая модель (учет гравитационных возмущений от больших планет, Земли и Луны, трех малых планет и релятивистских возмущений).
Таблица 7.2. Минимальные расстояния между Землей и Апофисом 13 апреля 2029 г.
В случае решения ИПА добавлен учет светового давления и сжатия Земли и Солнца, а также учет эффекта фазы. Модели движения, использованные в двух других случаях, лишь незначительно отличались от модели, использованной в ИПА. Несмотря на некоторое различие использованных моделей, видно, что все четыре решения дают близкие значения минимального расстояния между Землей и Апофисом.
Следующим этапом рассмотрения сближения является построение эллипса рассеяния в плоскости цели (см. раздел 7.5). Расчеты, выполненные на момент времени JD 2462240,407115, дают ξ = 7125 км, ζ = 37 550 км, aξ = 15,0 км, aζ = 351,6 км (рис. 7.5).
Если эллипс рассеяния ограничить только теми точками, которые образованы виртуальными астероидами, берущими начало в области начального эллипсоида рассеяния с осями ±3σEi, то его большая полуось равна 3aζ, а малая — 3σEi. Эллипс рассеяния, таким образом, весьма вытянут (рис. 7.5). Расстояние от центра эллипса рассеяния до центра Земли столь велико и размеры эллипса столь малы, что ни о каком столкновении в 2029 г. речи быть не может.
Рис. 7.5. Положение эллипса рассеяния на плоскости цели 13 апреля 2029 г.
7.7.4. Возможность резонансных возвращений в 2036 г. и в последующий период. Различные точки большой полуоси эллипса рассеяния соответствуют различным виртуальным астероидам. Если в номинальном решении варьировать значение среднего движения в пределах от +3σn до величины -3σn, а все остальные элементы оставлять неизменными, то точки пересечения виртуальных астероидов с плоскостью цели пробегают всю большую ось эллипса, начиная с ближайшего к Земле ее конца и заканчивая наиболее удаленным концом (рис. 7.5). Соответствующие виртуальные астероиды пройдут на различных расстояниях от центра Земли и поэтому их орбиты изменятся по-разному.
Особенно значительны по своей величине и возможным последствиям будут изменения большой полуоси. Так, среднее движение астероида, прошедшего через ближайший к Земле конец большой оси эллипса, изменится от 1,11385 до 0,84407 °/сут, а среднее движение астероида, прошедшего через дальний конец большой оси эллипса, изменится от 1,11385 до 0,85429 °/сут. Измененным значениям среднего движения соответствуют периоды обращения, выраженные в годах, P = 1,1677 и P = 1,1537. В силу непрерывности существуют виртуальные астероиды, которые будут иметь периоды обращения, равные любому числу в указанных пределах. В частности, верхний предел близок к отношению 7: 6 ≈ 1,1667. Виртуальные астероиды с периодами, близкими к 1,1667 года, по истечении семи лет, совершив шесть оборотов вокруг Солнца, опять окажутся вблизи Земли, и минимальное расстояние от Земли для некоторого множества из них может оказаться меньше или равным радиусу Земли, что будет означать столкновение. Более точное представление о реальной ситуации в апреле 2036 г. можно получить, если численным путем проследить движение большого числа виртуальных астероидов, чьи точки пересечения с плоскостью цели в апреле 2029 г. располагаются вдоль большой оси эллипса рассеяния. Получить начальные условия для таких виртуальных астероидов можно путем варьирования среднего движения Апофиса в начальную эпоху в пределах ±3σn (для большей гарантии вариацию можно брать в более широких пределах). Значения минимальных расстояний между Апофисом и Землей в 2036 г., полученные для вариаций среднего движения в пределах от +11 252 10-11 до +11 257 10-11 °/сут, приведены в табл. 7.3.
Таким образом, вариация номинального значения среднего движения примерно в указанных пределах ведет к столкновению Апофиса с Землей в 2036 г. Так как согласно решению ИПА σn = 2739 10-11 °/сут, то указанные вариации, хотя они и несколько превышают величину 4σn, тем не менее, имеют не исчезающе малую вероятность. Интересно выяснить, в каком месте плоскости цели соответствующие виртуальные астероиды пересекают ее. Очевидно, что точки пересечения располагаются на продолжении большой оси эллипса рассеяния, занимая некоторый его отрезок. Один конец отрезка соответствует величине вариации 4,1082σn, а второй — вариации, равной 4,1100 σn. На плоскости цели вариация, равная σn, соответствует величине aζ = 351,6 км. Из этого следует, что один конец отрезка находится на расстоянии 4,1082 351,6 = 1444,4 км от центра эллипса, а другой — на расстоянии 4,1100 351,6 = 1445,0 км, т. е. длина отрезка составляет около 600 м (рис. 7.5).
Таблица 7.3. Минимальные расстояния между Апофисом и Землей в 2036 г. при различных вариациях среднего движения астероида (решение ИПА)
Это так называемая «замочная скважина» (англ. keyhole) на плоскости цели, через которую должен пройти центр Апофиса, чтобы через шесть лет астероид столкнулся с Землей. Вариация других элементов орбиты астероида мало влияет на этот вывод. С учетом их вариации отрезок большой оси превратится в щель на плоскости цели, ограниченную с боков контуром эллипса шириной около 600 м, но результат в принципе останется прежним.
Совершенно аналогично для решения [Giorgini et al., 2008] при указанных в табл. 7.4 значениях вариации среднего движения находим соответствующие им значения Δmin (номинальное значение среднего движения равно 1,1128077236422 ± 1733 10-11 °/сут, эпоха JD 2453979,5).
Таблица 7.4. Минимальные расстояния между Апофисом и Землей в 2036 г. при различных вариациях среднего движения астероида (решение [Giorgini et al., 2008])
Итак, уже при иных значениях вариации среднего движения решение [Giorgini et al., 2008] также приводит к столкновению, хотя и в чуть более поздний момент времени. В табл. 7.5 сопоставлены расчеты вероятности столкновения Апофиса с Землей в 2036 г., выполненные на основе различных решений.
В таблице для каждого из трех решений последовательно приводятся следующие данные: cреднее движение в исходную эпоху и оценка его ошибки; минимальное расстояние астероида от Земли в апреле 2029 г.; интервал вариаций среднего движения, приводящих к столкновению в апреле 2036 г.; варьированное значение среднего движения в эпоху JD 2454200,5, соответствующее средней по величине вариации; минимальное расстояние от центра Земли в 2036 г.; расстояние от центра Земли до середины «замочной скважины» (скважины для столкновения, или зоны резонансного возврата), ведущей к столкновению в 2036 г.; ширина скважины, выраженная в единицах ошибки среднего движения для данного решения; вероятность столкновения Апофиса с Землей в 2036 г.
Таблица 7.5. Вероятность столкновения Апофиса с Землей согласно различным решениям
Следует особенно подчеркнуть два момента. Первый момент — это близкие положения «замочных скважин» для всех трех решений: разница их расстояний от Земли находится в пределах 2–3 км.
Второй момент, на который следует обратить внимание, — это очень большой разброс в оценке вероятности столкновения, вытекающей из различных решений. Оценка вероятности столкновения с учетом неполноты наших знаний о величине сил, действующих на астероид, будет рассмотрена в следующем параграфе.
Если виртуальный астероид пройдет мимо «замочной скважины», но близко к ней, то в 2036 г. столкновение будет исключено, но в результате такого тесного сближения с Землей орбита астероида претерпит новые изменения, которые могут привести к опасным сближениям в обозримом будущем. Это весьма неприятная ситуация, так как в результате двух последовательных тесных сближений астероида с Землей точность предсказания его дальнейшего движения резко снижается. Самые, на первый взгляд, малозначительные изменения действующих на астероид сил или потеря точности при интегрировании могут вести к неправильному прогнозу его движения после серии последовательных тесных сближений.
Скважины на плоскости цели 2029 г. отнюдь не единственные в своем роде (см. рис. 10.13 в главе 10). Например, на расстоянии всего в 142 км от центра эллипса рассеяния, по другую сторону от «замочной скважины», ведущей к столкновению в 2036 г., находится скважина для столкновения с Землей в 2051 г. Правда, в данном случае расчет показывает не прямое столкновение, а лишь касательное прохождение астероида мимо Земли на высоте всего в 2000 км. Ширина «замочной скважины», ведущей к касательным прохождениям на высотах, не превышающих 3500 км над поверхностью Земли, составляет 165 м.
7.7.5. Иcследование резонансных возвратов с помощью метода точечных гравитационных сфер. Как уже обсуждалось в этой главе, каждое тесное сближение астероида с планетой сопровождается рассеянием трубки его возможных траекторий и заметной потерей точности прогнозирования движения астероида. Траектории виртуальных астероидов, близкие до прохождения около планеты, претерпевают различные возмущения и расходятся после сближения с ней. С аналогичной трудностью небесные механики сталкивались уже давно при исследовании комет, сближающихся с Юпитером. Она имеет место и для астероида Апофис, а также для других опасных астероидов в случае их тесных сближений с планетами.
После сближения Апофиса с Землей в 2029 г. точность прогнозирования его движения резко падает. Если же для части возможных траекторий после 2029 г. будет иметь место тесное сближение (не соударение!) в 2036 г., то дальнейшее движение по этим траекториям становится практически недетерминированным. Другими словами, современной точности начальных данных недостаточно для исследования методом численного интегрирования возможных соударений в период, последующий за 2036 г., если в этом году будет иметь место тесное сближение астероида с Землей. По мере получения новых наблюдений точность орбиты будет, конечно, расти, однако в ближайшие несколько лет вряд ли удастся достаточно точно предвычислить траектории тесных сближений в 2036 г., чтобы надежно проанализировать все возможные сближения и возможные столкновения в последующие годы.
Интуитивно ясно, что опасными траекториями, на которые может перейти астероид в результате тесного сближения с Землей, являются резонансные траектории невысокого порядка. Например, если после сближения астероид окажется на орбите, которая имеет соизмеримость средних движений 1:1 с Землей, то ровно через год астероид и Земля окажутся в том же месте и возможно их новое тесное сближение или соударение. Из-за расширения трубки возможных траекторий в результате взаимодействия с Землей существует больше возможностей попадания в такие опасные резонансы. Соответствующие области начальных условий называют зонами резонансного возврата. Так, в работе [Chesley, 2005] рассмотрены возможные резонансные возвраты астероида Апофис после сближения в апреле 2029 г. Ко времени написания работы Чизли орбита астероида имела еще сравнительно низкую точность, и в результате рассеяния трубки возможных траекторий после этого сближения астероид мог оказаться на различных резонансных орбитах, указанных на рис. 10.13 в главе 10. Эти орбиты приводили к опасным сближениям Апофиса с Землей в 2034, 2035, 2036, 2037, 2038, 2046 и 2048 гг. Однако по мере увеличения точности орбиты число возможных резонансных возвратов быстро уменьшалось, и уже к августу 2005 г. остался лишь один возможный резонансный возврат в апреле 2036 г. Существенно в этой связи, что в ближайшие несколько лет до 2012 г. точность орбиты астероида Апофис не будет расти в результате поступления новых наблюдений, поскольку в это время наблюдать его с Земли практически невозможно. После появления новых наблюдений и уточнения орбиты Апофиса в 2012–2013 гг. совершенно аналогично уменьшится число возможных траекторий с опасными сближениями в период после 2036 г., которые детально обсуждаются в работе [Соколов и др., 2008].
В первом приближении резонансные возвраты можно исследовать, используя прием, называемый методом точечных гравитационных сфер (ТГС). Для нескольких последовательных сближений с планетой этот метод позволяет с умеренной трудоемкостью получить картину возможных преобразований орбит. При этом, если сближения достаточно тесные и относительная скорость достаточно велика, метод дает приемлемую точность результатов. Точность определяется из сравнения с результатами численного интегрирования неупрощенных уравнений движения. Метод ТГС позволяет в ряде важных случаев рассматривать вместо сложной задачи трех тел последовательность простых задач двух тел. Этот прием в различных модификациях часто применяется в небесной механике и механике космического полета, в частности для проектирования траекторий космических аппаратов с несколькими гравитационными маневрами. Некоторые подробности, а также литературные ссылки можно найти, например, в работах [Соколов и др., 2008; Елькин и др., 2003].
В методе ТГС исследуемая траектория аппроксимируется последовательностью кеплеровых гелиоцентрических орбит соударения с планетами. При «соударении» происходит мгновенное преобразование скорости — поворот вектора планетоцентрической скорости на угол между асимптотами планетоцентрической гиперболы. Модуль планетоцентрической скорости остается при этом постоянным. Элементы кеплеровой орбиты преобразуются в этом случае по формулам классической задачи двух тел. Как известно, по координатам и скоростям в данный момент времени элементы орбиты определяются однозначно. Угол поворота вектора планетоцентрической скорости зависит от прицельного расстояния и в рамках рассматриваемой аппроксимации (сжатия гравитационной сферы в точку, обнуления прицельного расстояния) не определен, т. е. может быть произвольным в некоторых пределах.
Для исследования траекторий, содержащих опасные сближения и соударения астероида Апофис с Землей сразу после возможного тесного сближения в 2036 г., в работе [Соколов и др., 2008] используются следующие методы. В предположении, что движение в промежутке между 2029 и 2036 гг. происходит по резонансной орбите соударения с соизмеримостью средних движений 6:7, методом ТГС строится множество возможных резонансных кеплеровых орбит соударения после тесного сближения в 2036 г. Затем с использованием численного интегрирования уравнений движения с учетом всех возмущений ищутся траектории тесных сближений и соударений с Землей, близкие к построенным кеплеровым резонансным орбитам соударения. Для преодоления основной трудности — потери точности из-за тесных сближений, область допустимых начальных данных можно транспортировать вдоль траекторий в 2035 г., в результате чего ее размеры увеличиваются на несколько порядков. После этого численное нахождение начальных данных, соответствующих опасным траекториям, не представляет принципиальных трудностей. Полученные методом ТГС приближенные значения минимальных геоцентрических расстояний используются при нахождении интересующих нас опасных сближений в качестве первых приближений, начальные данные для численного интегрирования берутся из новой области в 2035 г. Аналогичный подход использовал Ж. Ласкар при исследовании динамического хаоса в Солнечной системе и поиске «уходящих» из нее траекторий Меркурия.
В результате с использованием методов, указанных в работе [Соколов и др., 2008], были построены порождающие эллиптические резонансные траектории соударений в 2037–2052 гг. (95 траекторий). Некоторые из них не могут быть реализованы, так как проходят в 2036 г. на геоцентрическом расстоянии менее радиуса Земли, однако большинство соответствует реальным траекториям. Подробнее были рассмотрены порождающие траектории с соударениями до 2041 г. включительно: численно получены траектории, содержащие соударения Апофиса с Землей в 2040 и 2041 гг. (резонансы 3:4 и 6:5 после 2036 г.), а также тесные сближения в 2037, 2038 и 2039 гг. Численное интегрирование подтвердило корректность использования метода ТГС для построения приближенных решений в случае Апофиса. Для численного интегрирования использовался интегратор Эверхарта. Движение планет и Луны описывалось известными динамическими моделями DE403 и DE405. Минимальные геоцентрические расстояния после 2036 г., полученные по DE403 и DE405, различаются незначительно, разница обычно менее 1 км. Области начальных значений координат в 2035 г., соответствующие указанным соударениям в 2040 и 2041 гг., имеют размеры порядка сотен и десятков метров соответственно. Отметим, что размеры областей начальных данных, соответствующие соударениям, по DE403 и DE405 практически совпадают, в то время как сами эти области отстоят друг от друга на десятки километров.
Представляет также интерес история сближений Земли и Апофиса. Проведенное численное интегрирование от 2006 г. до 1700 г. показало отсутствие тесных сближений на этом интервале. Все сближения происходят 12–14 апреля. Самое тесное сближение имело место в 1819 г. до расстояния 0,84 млн км. Остальные шесть сближений происходили до расстояний 3–4 млн км.
На примере астероида Апофис хорошо видны некоторые важные аспекты проблемы астероидно-кометной опасности, на которые раньше не обращали должного внимания. Так, соударению астероида с планетой могут предшествовать сближения с ней, хотя бы потому, что сближения на умеренные расстояния более вероятны, чем соударения. Этот вопрос обсуждается, в частности, в работе [Елькин, Соколов, 1995]. Отмеченное обстоятельство позволяет заблаговременно обнаружить потенциально опасный объект. С другой стороны, рассеяние возможных траекторий при тесных сближениях ведет к потере точности, затрудняет прогнозирование движения и требует применения специальных методов. Среди возможных движений астероида после тесных сближений появляются опасные траектории, имеющие резонансные возвраты к Земле. Соответствующие «замочные скважины», или зоны резонансного возврата, имеют очень малые размеры и, следовательно, мала вероятность столкновения с Землей. Точное положение «скважин» зависит от величины возмущающих сил: малые недостаточно известные эффекты могут изменить место «скважин» в пространстве начальных данных. Точное определение их положения является сложной актуальной задачей. В то же время само наличие «скважин» и их размеры мало зависят от возмущений. Следует отметить сложную структуру соответствующего множества, аналогичную фрактальной. Небольшое искусственное изменение траектории астероида (удар по нему и т. п.) с целью предотвратить соударение в апреле 2036 г. не гарантирует отсутствия соударений в последующие несколько лет. Отклонив траекторию, можно попасть в близлежащую «замочную скважину», ведущую к соударению. Недостаточная точность знания орбиты Апофиса не позволяет пока исключить возможность попадания в эти зоны резонансного возврата. Поэтому наряду с совершенствованием методов прогнозирования движения астероидов особую ценность имели бы наблюдения АСЗ из космоса, а также использование для уточнения орбиты Апофиса сигналов радиопередатчика, доставленного на орбиту искусственного спутника этого астероида.
7.7.6. Влияние эффекта Ярковского на движение Апофиса. Во всех рассмотренных случаях модель движения Апофиса не включала влияние эффекта Ярковского (см. раздел 3.6). Он вызывается неравномерным нагревом тела солнечными лучами в результате осевого вращения тела и его движения по орбите. Эффект Ярковского зависит от положения оси вращения, орбиты и массы тела, от теплопроводности его поверхностных слоев. Поскольку большинство этих параметров неизвестны, явным образом учесть эффект невозможно. Однако неявным образом этот эффект частично учитывается в результате подгонки элементов орбиты, прежде всего большой полуоси, к наблюдениям, которые, естественно, отражают влияние эффекта, если он достаточно велик. Тем не менее, значительная часть эффекта остается неучтенной, что может вызвать заметные ошибки при прогнозировании движения тела на основе найденной орбиты. В работе [Giorgini et al., 2008] выполнена оценка максимального смещения Апофиса по орбите за период от эпохи оскуляции (2006 г.) до 2029 г. под действием не учитываемой части эффекта. При предполагаемой массе в зависимости от направления вращения и гипотетически заданной теплопроводности поверхностных слоев максимальное смещение вдоль орбиты может составить к 2029 г. от –720 км до +780 км. Эти значения находятся в качественном согласии с оценками влияния эффекта Ярковского на движение ряда АСЗ, найденными Ю. А. Чернетенко иным путем по сравнению с работой [Giorgini et al., 2008].
В работе [Чернетенко, 2007] учет эффекта Ярковского в движении астероидов производился без каких-либо предположений о физических характеристиках этих тел. Предполагалось лишь, что зависимость этого эффекта от гелиоцентрического расстояния r имеет вид 1/r2, а его величина характеризуется в общем случае тремя составляющими ускорения: радиальной, трансверсальной и нормальной, которые определяются из наблюдений (радарных и оптических) совместно с параметрами орбиты. При этом для астероида (6489) Голевка включение в число определяемых по наблюдениям запаздывания параметров трансверсальной составляющей A2 позволило уменьшить ошибку единицы веса с 2,3 мкс до 0,5 мкс. Величина A2 оказалась равной (-2,00 ± 0,14) 10-14 а.е./сут2.
При оценке влияния эффекта Ярковского на движение астероида Апофис было принято, что зависимость ускорения от гелиоцентрического расстояния имеет вид 1/r2, а для A2 принимались некоторые возможные значения (+2, –2, +6, –6) 10-14 а.е./сут2. Значения минимального расстояния от Апофиса до Земли и смещения астероида вдоль орбиты на 13 апреля 2029 г., полученные в результате уточнения параметров движения и последующего интегрирования, приведены в табл. 7.6.
Таблица 7.6. Влияние принятых значений ускорения, вызываемого эффектом Ярковского, на величину минимального расстояния от Апофиса до Земли и смещение астероида вдоль орбиты на 13 апреля 2029 г.
Изменения смещения вдоль орбиты при A2 = (+6,-6) 10-14 а.е./сут2 близки к оценкам, полученным в работе [Giorgini et al., 2008]. Конечно, подобные расчеты имеют только качественный характер, поскольку реальное значение A2 из имеющихся наблюдений найти невозможно. Но метод может быть использован в будущем, если будет получено необходимое число достаточно точных наблюдений.
Оценим теперь, как влияет эффект Ярковского на положение эллипса рассеяния на плоскости цели и на вероятность столкновения Апофиса с Землей в 2036 г. Будем исходить из предположения, что величина A2 для астероида равна -2 10-14 а.е./сут2, что приводит к сокращению его полуоси и к увеличению среднего движения. В результате к апрелю 2029 г. астероид сместится вдоль орбиты на величину +272 км (табл. 7.6). Таким образом, астероид, соответствующий номинальному решению, найденному с учетом дополнительного ускорения, пересечет плоскость цели несколько раньше, и его минимальное расстояние от Земли составит, как показано в табл. 7.6, 37 964 км по сравнению с расстоянием в 38 220,5 км, на котором астероид должен проследовать мимо Земли в 2029 г. согласно решению, полученному ИПА.
7.7.7. Перспективы уточнения орбиты Апофиса. Для планирования возможного изменения орбиты Апофиса в случае, если это изменение будет вести к неизбежному или весьма вероятному столкновению, важно иметь представление о том, каковы перспективы проведения дальнейших астрометрических и физических наблюдений данного астероида и какого уточнения модели движения астероида можно ожидать в предстоящий период. При рассмотрении данного вопроса важно иметь в виду, что активные действия по изменению орбиты Апофиса должны быть предприняты до 2029 г. Величина корректирующего одиночного импульса, прикладываемого к астероиду до апреля 2029 г., чтобы избегнуть столкновения в 2036 г., примерно на три порядка меньше необходимого импульса после сближения в 2029 г. [Ивашкин, Стихно, 2008]. Причина такого различия вполне понятна. В случае, если коррекция выполняется до 2029 г., необходимо лишь несколько изменить минимальное геоцентрическое расстояние астероида (изменить орбиту так, чтобы ее пересечение с плоскостью цели произошло вне ближайшей окрестности «замочной скважины» размером около 600 м). Остальную работу по изменению траектории выполнит гравитационное поле Земли. После 2029 г. гелиоцентрическую орбиту астероида придется смещать уже на величину порядка нескольких радиусов Земли.
Размеры астероида и особенности его орбиты таковы, что благоприятные условия для его наблюдения складываются лишь в определенные весьма непродолжительные периоды времени. На рисунках 7.6–7.8 показаны элонгация астероида Апофис от Солнца, его звездная величина и расстояние от Земли в период 2006–2024 гг. На протяжении 2009–2011 гг. видимая звездная величина астероида будет оставаться в пределах 19,5m–22,1m. Но относительно большой блеск имеет место при малых элонгациях от Солнца. Сомнительно, что при звездной величине, превышающей 21m, астероид можно будет наблюдать в элонгациях, меньших 50°. Таким образом, ближайшее «окно» для оптических наблюдений откроется не ранее декабря 2011 г., когда элонгация превысит 50°, а звездная величина будет меньше 21m.
Рис. 7.6. Элонгация астероида Апофис от Солнца
Возможности проведения новых радиолокационных наблюдений Апофиса также очень ограничены. Для наиболее крупного 300-м радиолокатора в Аресибо (Пуэрто-Рико, США) предельно возможное расстояние приема отраженного от астероида сигнала (при заданных размерах астероида и значении радиоальбедо его поверхности) составляет 0,31 а.е. Для 70-м локатора в Голдстоуне (США) расстояние приема составляет всего 0,14 а.е.
Рис. 7.7. Видимая звездная величина астероида Апофис
Рис. 7.8. Расстояние от астероида Апофис до Земли
Примерно такие же возможности имеются у радиолокатора в Евпатории и модернизируемого радиотелескопа в Уссурийске. Таким образом, «окна» для радиолокационных наблюдений будут открыты только в конце 2012 г. — первой половине 2013 г., в октябре 2020 г. и феврале — августе 2021 г.
В работе [Chesley, 2005] предпринята попытка оценить, насколько точно может быть определена орбита Апофиса из возможных оптических и радиолокационных наблюдений, при условии, что из них будет получена не только астрометрическая информация о положениях Апофиса, но и найдено положение оси и направление вращения астероида. Как отмечалось ранее, две последние величины важны для оценки влияния эффекта Ярковского на движение астероида. В работе предполагается, что в периоды, наиболее благоприятные для наблюдений, каждую вторую ночь будет получено одно астрометрическое наблюдение со средней ошибкой в 0,2″, а в периоды, когда наблюдения будут менее доступны, но все еще возможны для крупных телескопов, одно наблюдение с той же ошибкой будет получено в течение каждой лунации. Далее предполагается, что в периоды 14–20 февраля и 6–10 июля 2013 г., а затем в периоды 9–12 октября 2020 г. и 16–20 марта 2021 г. с помощью радиолокатора в Аресибо будет выполнено в сумме 23 измерения расстояния до астероида с различной точностью в зависимости от расстояния. С учетом уже имеющихся наблюдений и предположительно осуществленных к исходу того или иного интервала можно оценить, с какой точностью определяется большая полуось эллипса рассеяния на плоскости цели в 2029 г. по всем имеющимся на рассматриваемый момент наблюдениям. Из вычислений Чизли следует, что к началу 2014 г. большая полуось эллипса рассеяния сократится примерно до 30 км, а к началу 2022 г. — до величины, несколько превышающей 10 км, если радарные наблюдения осуществлены не будут, и до 1 км, если они будут выполнены. Кроме того, в работе Чизли моделируется уточнение эллипса рассеяния с учетом годичного приема радиосигналов с эквивалентной точностью ±2 м от передатчика, доставленного на Апофис в 2019 г. С учетом радионаблюдений точность определения большой полуоси эллипса рассеяния повысится до 100 м.
В работе [Giorgini et al., 2008] эта проблема также обсуждается, но с применением иных подходов. Впрочем, результаты в целом оказываются сопоставимыми.
Можно поставить вопрос: какова вероятность того, что с учетом всех факторов столкновение с Апофисом в 2036 г. произойдет, и как вероятность столкновения будет эволюционировать с течением времени? Можно попробовать найти ответ на этот вопрос, если еще раз обратиться к рис. 7.5. Мы уже видели, что «замочная скважина» для столкновения в 2036 г. располагается на расстоянии 1444,7 км = 4,11σζ от центра эллипса, что определяет малую вероятность столкновения. Но при этом не был принят во внимание эффект Ярковского. Он способен сдвинуть центр эллипса, а вместе с ним и весь эллипс рассеяния вдоль большой оси в ту или иную сторону. Направление смещения и его величина зависят от положения оси вращения астероида. Если наклон оси вращения к плоскости орбиты астероида меньше 90° (прямое вращение), то, как нетрудно видеть, реактивный эффект покидающих тело тепловых фотонов имеет составляющую, направленную по вектору орбитальной скорости тела (рис. 7.9).
Рис. 7.9. Эффект Ярковского в случае прямого вращения астероида; F — вектор реактивного ускорения, v — орбитальная скорость астероида
Это приводит к увеличению большой полуоси орбиты астероида и уменьшению его среднего движения. Астероид будет двигаться по орбите с некоторым отставанием во времени и догонит плоскость цели 13 апреля 2029 г., несколько позже по сравнению с невозмущенным случаем. Из этого следует, что координата ζ центра эллипса окажется больше и астероид пройдет мимо Земли на несколько большем расстоянии от ее центра и от «замочной скважины» (положение последней в возмущенном случае останется практически тем же). В случае обратного осевого вращения астероида все будет обстоять с точностью до наоборот. Астероид пересечет плоскость цели в более ранний момент времени по сравнению с невозмущенным случаем, и весь эллипс как целое сместится в сторону меньших геоцентрических расстояний. По расчету [Giorgini et al., 2008] эффект Ярковского может привести к смещению астероида вдоль орбиты на величину от 325 км до 740 км в зависимости от принятых значений радиуса, плотности вещества и коэффициента теплопроводности. Можно рассчитать, что значение координаты ζ при этих смещениях вдоль орбиты изменится на величины того же порядка. При прямом направлении вращения это означает, что расстояние от центра эллипса до «замочной скважины» возрастет на указанные величины и что столкновение практически невозможно (вероятность столкновения меньше 10-7). При обратном направлении вращения эффект Ярковского способен сдвинуть эллипс рассеяния так, что «замочная скважина» окажется в пределах 3aζ. Вероятность столкновения в этом случае остается малой, но заметно отличной от нуля.
7.8. Астероид 2008 ТС3
Данный объект называют астероидом, поскольку он был обнаружен по технологии наблюдения астероидов и получил первичное обозначение как астероид. На самом деле это тело оказалось малым по размеру и, строго говоря, должно называться не астероидом, а метеороидом (см. главу 5). Тем не менее, этот объект интересен тем, что его столкновение с Землей было заранее спрогнозировано и реально наблюдалось. Тем самым международная кооперация по обнаружению АСЗ прошла важное испытание в отношении оперативности и слаженности ее действий, а обсуждающиеся в этой главе методы прогноза столкновений прошли практическую проверку. Все события с момента открытия астероида и до момента его столкновения с Землей произошли очень быстро и заняли не более суток.
6 октября 2008 г. в 6 ч 39 мин Всемирного времени (UTC) Р. Ковальский на обсерватории Маунт Леммон (штат Аризона, США) обнаружил объект примерно 19m, предварительно обозначенный им как STA9D69. Наблюдения выполнялись в рамках финансируемого НАСА обзора объектов, сближающихся с Землей (Catalina Sky Survey). Наблюдения немедленно были переданы в Международный планетный центр, где была вычислена предварительная орбита, показавшая, что объект с вероятностью 98 % находится на орбите, приводящей к столкновению с Землей в ближайшее время. Объект получил предварительное обозначение 2008 ТС3. Планетный центр немедленно разместил информацию об объекте на своем сайте и передал ее в Лабораторию реактивного движения НАСА (США). В течение часа после получения первых данных в ЛРД было вычислено, что объект войдет в атмосферу Земли 7 октября над северо-восточной Африкой.
Исходя из оценок звездной величины и расстояния до астероида его размер был определен в 2–5 м. Обычно объекты такого размера полностью разрушаются в атмосфере Земли. Поверхности Земли достигают только небольшие осколки, потерявшие высокую скорость. Вхождение астероида в атмосферу могло наблюдаться в течение нескольких десятков секунд как яркий болид. По имеющимся данным, вторжения тел таких размеров в атмосферу Земли происходят несколько раз в год. Следовательно, он не представлял особой угрозы для Земли. Однако оповещенные наблюдатели поторопились к телескопам, чтобы наблюдать этот необычный астероид до того, как он войдет в тень Земли и станет невидимым. В этих наблюдениях приняли участие и российские наблюдатели: астроном-любитель Т. Крячко, наблюдавший астероид на Кавказской астрономической станции обсерватории им. В. П. Энгельгардта Казанского государственного университета. На автоматизированном телескопе ГАО РАН (Пулково) ЗА-320М в ночь с 6 на 7 октября 2008 г. были проведены астрометрические и фотометрические наблюдения астероида 2008 ТС3. За 4 часа получено 270 наблюдений в интегральной полосе телескопа. Оценена абсолютная звездная величина астероида, его размер (4,8 ± 0,8 м) и масса (131 ± 5 т).
Новые наблюдения позволили уточнить первоначальную орбиту астероида. В течение последующих 11 часов было выпущено 24 электронных циркуляра, содержащих уточняемые параметры орбиты и эфемериды астероида. Оказалось, что большая полуось орбиты составляет 1,27 а.е., эксцентриситет — 0,29, а период обращения — 1,43 года, q = 0,91 а.е., Q = 1,63 а.е. Стив Чизли (ЛРД) сообщил, что астероид войдет в земную атмосферу примерно в 2 ч 45 мин 28 с UTC и достигнет максимального торможения на высоте примерно в 14 км в 2 ч 45 мин 54 с UTC с неопределенностью в ±15 с. Всего на 27 обсерваториях было получено 589 позиционных наблюдений; блеск астероида, зафиксированный в последних наблюдениях, составил примерно 13–14m.
А. Харрис, М. Козубал, Р. Дантовиц и П. Правек из наблюдений и анализа кривой блеска определили, что астероид совершает сложное вращение вокруг двух осей с периодами 97 с и 49 с.
Астероид вошел в атмосферу со скоростью примерно 12,4 км/с. Энергия взрыва оценивается примерно в одну килотонну в тротиловом эквиваленте. Это произошло над малонаселенным районом Земли (северный Судан), и, возможно, единственными непосредственными наблюдателями этого события были заранее предупрежденные о нем члены экипажа самолета авиалинии Эйр Франс — KLM, который находился от предполагаемого места входа астероида в плотные слои атмосферы на расстоянии примерно в 1400 км. Пилоты наблюдали короткую яркую вспышку примерно в предсказанное время. Изображение высотного следа в атмосфере было получено на севере Судана на рассвете 7 октября (см. рис. 7.10 на вклейке).
Инфразвуковой детектор, расположенный в Кении, зарегистрировал звуковую волну, исходящую из предполагаемого района падения. Она соответствовала энергии 1,1–2,1 килотонны в тротиловом эквиваленте. Cпутник Метеосат-8, выполнявший наблюдения этого события в визуальном и инфракрасном диапазонах, зафиксировал короткую вспышку в атмосфере над районом предполагаемого падения астероида (рис. 7.11).
Эта короткая драматическая история астероида 2008 ТС3 показала, что существующие службы слежения за ОСЗ в состоянии обнаружить объект на опасной для Земли траектории на его подлете к Земле, а вычислительные центры — оперативно обработать наблюдательную информацию и информировать заинтересованных лиц о предстоящих событиях.
Рис. 7.11. Вспышка в атмосфере над северным Суданом, зафиксированная спутником Метеосат-8 над районом предвычисленного места падения астероида 2008 TC3 (http://www.eumetsat.int/Home/Main/News/Features/707785/). Copy-right @ EUMETSAT (European Organisation for the Exploitation of Meteorological Satellites)
В течение нескольких часов после обнаружения объекта весь мир уже знал, где астероид можно наблюдать, где и когда он упадет, а также представляет ли это угрозу для землян или нет. Однако возникает вопрос: «А что бы произошло, если бы астероид был гораздо больше по размеру?» Можно только сказать, что более крупный объект может быть обнаружен не за 20 часов до столкновения, а раньше. Так, объект с диаметром в несколько десятков метров может быть обнаружен за несколько дней до столкновения. Тем не менее, этого времени явно недостаточно для предотвращения столкновения. Можно лишь предпринять меры по эвакуации населения из предполагаемого района падения.