Как помочь телескопу исследовать экзопланету?⁠⁠

В 2029 году Европейское космическое агентство запустит телескоп ARIEL для  исследования химического состава атмосфер экзопланет. Чтобы учесть возможные технические проблемы и решить их заранее, ежегодно проводится международное состязание Ariel Data Challenge на платформе Kaggle (онлайн-сообщество для проведения соревнований по Data Science). Студенты Института компьютерных наук Университета МИСИС Георгий Апарин и Артём Хапилов в этом году завоевали золото, решив кейс по улучшению качества данных, которые получают со спектрометров спутника. Георгий подробно рассказал о техническом решении задачи.
Гуманитарии, простите.

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

Это позволяет «поймать» фотоны, которые пролетели через атмосферу экзопланеты. Они взаимодействуют с химическими компонентами, теряя энергию в процессе. Разные молекулы поглощают свет на определенных длинах волн, что создает характерные особенности в спектре. Сравнивая фоновое излучение звезды до транзита и во время него, можно выявить особенности и определить, какие молекулы присутствуют в атмосфере.

На первый взгляд может показаться, что Y должна быть постоянной, ведь радиусы звезды и планеты неизменны. Есть небольшое «но». Радиус планеты в этом контексте зависит от длины волны. Определенные молекулы в атмосфере поглощают больше света, создавая эффект, будто перед звездой проходит планета с чуть большим радиусом, но без атмосферы. Это одна из особенностей астрономической спектроскопии, позволяющая изучать состав атмосфер, которые находятся от нас на огромном расстоянии.

На спутнике используются высокочувствительные спектрометры, настроенные на определенные длины волн, но работают они с погрешностью. Через атмосферу планеты пролетает от 50 до 200 ppm (фотонов на миллион) – крайне небольшое количество, которого недостаточно для прямой конвертации сырых данных в распределение. На Ariel будет установлено несколько оптических инструментов, специализирующихся на разных спектральных диапазонах и режимах наблюдения. FGS1 – это спектрометр системы точного наведения, задача которого обеспечить центрирование и фокусировку спутника, высокоточную фотометрию целевой звезды в видимом спектре. Он имеет чувствительность от 0,60 до 0,80 мкм. AIRS-CH0 – это инфракрасный спектрометр с чувствительностью от 1,95 до 3,90 мкм.

Задача соревнования – восстановить распределение из сырых данных спектрометров спутника при помощи различных методов. Все усложнялось наличием значительных шумов, среди которых особенно выделяется так называемый «джиттерный шум» (jitter noise) – результат микровибраций космического аппарата. Он создается вращением маховиков для стабилизации телескопа и схож с размытостью при попытке сделать фото с длинной выдержкой дрожащей камерой. Амплитуда джиттерного шума может достигать 200 ppm, что сопоставимо с величиной самого сигнала, обусловленного атмосферой планеты, и особенно затрудняет анализ сигналов от небольших планет.

Данные предоставлялись в виде временной последовательности 2D-изображений спектральной фокальной плоскости, симулируя реальные наблюдения телескопа. Необходимо было минимизировать шумы и извлечь из зафиксированной энергии значимую информацию для восстановления распределения. Для каждого наблюдения дан 3D-массив с осями (T, S, wl), где T – временная ось, S – пространственная, wl – спектральная. Результаты оценивались с помощью функции логарифмического правдоподобия Гаусса (Gaussian Log-Likelihood, GLL). Она сравнивает предсказанные спектр μ_user и отклонение σ_user с реальными данными по формуле:

Финальная метрика — отношение разности предсказанного GLL и «наивного» GLL к разности «идеального» GLL и «наивного» GLL, где наивный GLL представлял собой GLL с предсказаниями, равными среднему значению таргета по датасету (да, всё сложно))) В «идеальном» логарифмическом правдоподобии сигма была взята как 1e-5. Если мы предскажем значения, которые будут лежать в границах 1e-5 от таргета, метрика будет максимальной.

Стоит отметить, что на трейн данных даны наблюдения с двух различных звёзд – планетарных систем. В тесте к известным прибавлялись еще две неизвестные системы. Также распределение таргета на тесте отличалось от распределения на трейне, поэтому ни один подход, основанный на 2D или 3D-CNN не реализовался.

В нашем решении мы упростили постановку задачи. Данные имеют три размерности: время, пространство и спектр. Избавимся от пространственной размерности, усреднив всю энергию в одно значение для каждой длины волны. Для каждого наблюдения получим поверхность в координатах «время» и «длина волны». Для каждой пары (t_i, wl_i) у нас есть усредненное значение энергии – flux. Шум измерений высокий, поэтому для наглядности представлены графики одного и того же наблюдения с применением сглаживания и без него.

Обратите внимание на левую картинку: она визуально понятнее правой, поскольку на ней видно поверхность, составленную из нормализованных временных рядов для каждой длины волны относительно времени. Момент «провала» энергии — это момент транзита.

Теперь еще упростим задачу. Рассмотрим временной ряд для конкретной длины волны, вместо анализа всех:

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

Значения «до» и «после» отражают полную энергию звезды, а во время транзита – полную энергию за вычетом заблокированной. Если прибавить к значению кривой во время транзита заблокированную энергию, то получится состояние как будто планета не пересекала звездный диск. Отношение этого добавленного значения к фоновому – искомый таргет. Таким образом, верно восстановив полную энергию звезды на временном промежутке, мы сможем получить таргет. Для восстановления полной энергии была предложена идея полиномиальной интерполяции: полную энергию звезды с достаточной точностью можно описать полиномом.

Остаётся вопрос: как оценить лучший полином для интерполяции? Решение простое. Давайте считать MAE между предсказанным полиномом и кривой полной энергией звезды (далее – кривой светимости), при поднятии момента транзита на заданный коэффициент. Полином с наименьшей ошибкой будет соответствовать лучшему коэффициенту поднятия, которое по предварительной оценке и будет являться таргетом. Не решено, как формально поставить задачу, чтобы можно было применить к ней различные методы оптимизации. Обозначим кривую светимости с транзитом – S, влияние транзита на светимость – D. Также зададим гладкую функцию D через коэффициенты. Фазы начала и конца транзита – p1 и p2, время длительности транзитного перехода – коэффициент t, искомый таргет – коэффициент d. Поделив данный сигнал S на функцию мультипликации D, получим исходную кривую светимости, как на предыдущих графиках.

Для нахождения лучших коэффициентов был использован метод оптимизации Nelder-Mead в библиотеке scipy. В поставленной задаче мало параметров, поэтому он быстро сходится в глобальный оптимум. Оптимизатор делает свои функции, но оптимум по полиному не является решением задачи, а только упрощает его. Если зафиксировать коэффициент d и присвоить ему истинное значение, посчитать полином и MAE по восстановленной кривой светимости, то ошибка будет больше, чем при обычной оптимизации. Из этого следует, что в реальности переформулированная задача описывает таргет с достаточной точностью, но глобальный её оптимум не равен искомому таргету.

Сигналы для отдельных длин волн шумные, поэтому было решено объединить соседние длины волн особым образом и считать полином по усредненному сигналу.

Ниже приведён код основного алгоритма сбора предсказаний

1. Вычисляем p1, p2, и t по производной сигнала

2. Вычисляем d_st по усредненному сигналу. Этот параметр представляет собой предсказание среднего по таргету. Здесь мы также оптимизируем параметры, полученные на первом этапе.

3. Далее итерируемся по всем длинам волн с окном 2*k_binn_wl и шагомk_binn_wl // 2.

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

5. Далее каждый сабсемпл разбивается на подгруппы меньшего размера, для которых мы оптимизируем таргет отдельной функцией минимизации по их усредненному сигналу.

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

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

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