10 апреля 2019 года человечеству показали оранжевый бублик. Журналисты назвали его «первой фотографией черной дыры». Через час картинка была у всех — мемы про глаз Саурона, шутки про пончик, антропоморфизация,  заголовки «ученые сфотографировали невидимое».

Проблема в том, что это не совсем фотография.Точнее сказать, это очень странная фотография: если бы вы использовали телескоп горизонта событий (англ. EHT — далее по тексту) «как камеру» и нажали кнопку, вы бы получили черный квадрат и никакого бублика. Потому что он делает измерения, из которых алгоритм уже собирает изображение…  которого нет.

Вот про этот алгоритм и про то, как 3,5 петабайта данных летели в Бостон самолетом, и пойдет речь.

Зачем восемь телескопов, если можно один большой

Чтобы увидеть тень черной дыры в галактике M87 — той самой, в 55 млн световых лет от нас, — нужно угловое разрешение порядка 20 микросекунд дуги. Это, грубо говоря, как из Москвы разглядеть апельсин в Нью-Йорке с точностью до косточки.

Угловое разрешение телескопа определяется простой формулой:

Здесь λ — длина волны наблюдения, D — диаметр апертуры. EHT работает на 230 ГГц, то есть λ ≈ 1,3 мм. 

Подставим целевое разрешение:

13 000 км — это диаметр Земли (ну, чуть меньше на самом деле). Тарелку диаметром с Землю никто не построит, но есть другой путь.

VLBI: тарелка-призрак

Идея интерферометрии со сверхдлинной базой звучит так: если у нас есть два телескопа, разнесенные на расстояние B, то они вместе работают как кусок виртуальной тарелки диаметром B. Но не вся тарелка целиком — только одна, скажем так, «полоска» от нее.

Восемь телескопов EHT (от Чили до Южного полюса, от Мексики до Испании) дают C (8,2) = 28 парных базовых линий. Самая длинная — от Южного полюса до Испании — около 11 000 километров. Это и есть наш виртуальный диаметр. 

Но «28 полосок от тарелки» — это совсем не то же самое, что цельная тарелка. И тут пригождается наша любимая математика.

Облачная инфраструктура для ваших проектов

Виртуальные машины в Москве, Санкт-Петербурге и Новосибирске с оплатой по потреблению.

Подробнее →

Каждая пара телескопов измеряет одну точку Фурье-плоскости

Это ключевой момент, который ломает привычную логику (по крайней мере мою). Радиоинтерферометр измеряет Фурье-компоненты распределения яркости.

Если I (l, m) — это яркость неба в небесных координатах (l, m), то пара телескопов с базовой линией, проектируемой на небо как (u, v) (в единицах длины волны), измеряет величину, которая называется видимостью:

Это в точности двумерное Фурье-преобразование I. Каждая пара телескопов в каждый момент времени дает одну комплексную точку V(u, v).

Когда Земля вращается, проекция базовой линии на плоскость, перпендикулярную направлению на источник, меняется — каждая пара рисует на (u, v)-плоскости эллиптическую дугу. Это называется апертурный синтез: мы используем вращение Земли, чтобы заметать больше точек Фурье-плоскости, не двигая телескопы физически. 

99,9% данных отсутствует

Чтобы восстановить изображение по обратному Фурье, в идеале нужно знать V(u, v) во всех точках плоскости. У EHT же несколько изогнутых линий, нарисованных на огромной плоскости.

Если оцифровать (u, v)-плоскость в сетку 1 024×1 024, то у нас будет порядка 10⁶ ячеек. Реальные данные покрывают, дай бог, тысячу из них. Это и есть тот самый «99,9% данных нет» из заголовка.

Картинку покрытия легко получить руками:

import numpy as np
import matplotlib.pyplot as plt

plt.style.use('dark_background')

# Положения станций EHT (долгота, широта в градусах)
stations = {
    'ALMA': (-67.75, -23.03), 'APEX': (-67.76, -23.01),
    'JCMT': (-155.47, 19.82), 'SMA': (-155.48, 19.82),
    'SMT': (-109.89, 32.70), 'LMT': (-97.31, 18.99),
    'IRAM': (-3.39, 37.07), 'SPT': (45.00, -90.0),
}

R_earth = 6.371e6

def ecef(lon, lat):
    lon, lat = np.radians(lon), np.radians(lat)
    return R_earth * np.array([
        np.cos(lat) * np.cos(lon),
        np.cos(lat) * np.sin(lon),
        np.sin(lat)
    ])

coords = np.array([ecef(*ll) for ll in stations.values()])

# Все 28 базовых линий
n_stations = len(coords)
baselines = np.array([
    coords[i] - coords[j]
    for i in range(n_stations) for j in range(i + 1, n_stations)
])

# Эмулируем 6 часов наблюдения: вращаем Землю относительно источника
hours = np.linspace(-3, 3, 400) * np.pi / 12  # часовой угол в радианах
uv = []

for h in hours:
    Rz = np.array([
        [np.cos(h),  np.sin(h), 0],
        [-np.sin(h), np.cos(h), 0],
        [0, 0, 1]
    ])
    # Проекция базовых линий на плоскость UV
    uv.append((baselines @ Rz.T)[:, :2])

uv = np.concatenate(uv) / 1.3e-3  # переводим в длины волн при λ = 1.3 мм

# Визуализация
plt.figure(figsize=(7, 7))
plt.scatter(uv[:, 0], uv[:, 1], s=1, c='cyan')
plt.scatter(-uv[:, 0], -uv[:, 1], s=1, c='cyan')  # Эрмитова симметрия
plt.gca().set_aspect('equal')
plt.title('UV-coverage EHT (упрощенная модель)')
plt.xlabel('u, λ')
plt.ylabel('v, λ')
plt.grid(alpha=0.2)
plt.show()

Запустите и увидите редкие изогнутые лепестки на огромной черной плоскости. Вот это и есть «фото», которое мы получаем напрямую: где телескопы намерили, там данные есть. А между лепестками миллионы неизвестных Фурье-компонент.

Грязное изображение и грязный пучок

Что будет, если просто взять обратное Фурье от того, что есть? Формально наш набор измерений можно записать так:

Здесь S(u, v) — выборочная функция (sampling function), равная 1 там, где мы намерили, и 0 —там, где нет.

По теореме о свертке, обратное Фурье от произведения — это свертка обратных Фурье сомножителей:

Здесь B_dirty — обратное Фурье от S, известное как грязный пучок или PSF данной конфигурации. Это та функция размытия, через которую инструмент «смотрит» на небо.

И вот тут начинаются трудности: B_dirty выглядит ужасно. У него есть центральный пик, но вокруг — длинный шлейф боковых лепестков, артефактов и звезда Давида от регулярной структуры базовых линий. Если просто взять обратное быстрое преобразование Фурье, получится не M87 с тенью, а размазанная чертовщина, которую невозможно интерпретировать.

Нужна деконволюция. Нужно решить уравнение I_dirty = I_true * B_dirty относительно I_true. Только это уравнение не решается — оно недоопределено. Есть бесконечно много I_true, дающих то же самое I_dirty, просто потому что S(u, v) — нули в большинстве точек, и любая функция, которая в этих точках Фурье-плоскости делает что угодно, а в остальных совпадает с истинной видимостью, будет валидным решением.

Это классическая обратная задача и решается она, как все подобные задачи: введением априорных предположений о том, как должно выглядеть «правильное» изображение.

CLEAN: алгоритм 1974 года, который все еще работает

В 1974 году Ян Хегбом из Гронингенского университета предложил гениально простой алгоритм. Он называется CLEAN, и его до сих пор используют — в том числе одна из команд EHT.

Идея в следующем: предположим, что небо состоит из множества точечных источников. Если это так, то можно итеративно «выковыривать» их из грязного изображения.

import numpy as np
from scipy.signal import convolve

def clean(dirty_image, dirty_beam, gain=0.1, threshold=1e-3, max_iter=10000):
    residual = dirty_image.copy()
    components = np.zeros_like(dirty_image)
    
    # Центр грязного пучка для расчета сдвига
    center = np.array(dirty_image.shape) // 2

    for i in range(max_iter):
        # 1. Найти самый яркий пиксель в остатке
        peak_pos = np.unravel_index(np.argmax(np.abs(residual)), residual.shape)
        peak_val = residual[peak_pos]

        # 2. Если он сравним с шумом — стоп
        if np.abs(peak_val) < threshold:
            break

        # 3. Записать долю этого пика как «найденный точечный источник»
        components[peak_pos] += gain * peak_val

        # 4. Вычесть из остатка эту долю, свернутую с грязным пучком
        # Сдвигаем пучок так, чтобы его пик совпал с peak_pos
        shift = np.array(peak_pos) - center
        shifted_beam = np.roll(dirty_beam, shift, axis=(0, 1))
        
        residual -= gain * peak_val * shifted_beam

    # 5. Свернуть компоненты с «чистым пучком» (обычно Гауссиан) и добавить остаток
    # clean_beam должен быть определен заранее
    return components, residual

Что здесь происходит, по сути:

1. Мы предполагаем, что в самом ярком пикселе действительно сидит точечный источник;

2. Также, мы знаем, как один точечный источник выглядит после свертки с PSF — он выглядит как сам PSF, сдвинутый в эту точку;

3. Мы вычитаем небольшую долю этого «теоретического вклада» из изображения;

4. Повторяем тысячи раз;

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

CLEAN — алгоритм жадный, он работает потому, что для типичных радиоастрономических картинок (несколько ярких квазаров на пустом небе) предположение «все состоит из дельта-функций» оказывается отличной аппроксимацией.

Для тени черной дыры — уже совсем не так. Поверхность яркости вокруг M87* гладкая, ассимметричная, с кольцевой структурой. Хотя CLEAN все равно работает, но требует вариантов и тонкой настройки. Поэтому EHT использует и кое-что поновее.

Альтернатива — байесовский подход

Мы хотим найти такое изображение I, которое согласуется с измерениями и удовлетворяет нашим представлениям о том, что такое «разумная картинка».

Формально, минимизируем:

Первое слагаемое — это правдоподобие: насколько изображение согласуется с измеренными видимостями.

Здесь V_model — Фурье-преобразование текущего изображения I, посчитанное в тех точках, где мы намерили.

Второе слагаемое — это регуляризаторы. Их обычно несколько и они кодируют априорные предположения. Самые ходовые в EHT:

• гладкость: R_TV(I) = Σ |∇I| — штрафует резкие переходы;

• разреженность: R_ℓ₁(I) = Σ |Iᵢ| — штрафует ненужные ненулевые пиксели;

• энтропия (MEM): R_MEM(I) = −Σ Iᵢ · log(Iᵢ / mᵢ) — заставляет изображение быть «как можно более похожим на дефолтную модель m при равных χ²»;

• положительность: I ≥ 0 — отрицательной яркости не бывает;

• коэффициенты λ_k — гиперпараметры, которые нужно настраивать.

Слепые команды

У вас есть набор измерений с базовых линий. У вас есть алгоритм с десятком гиперпараметров. У вас есть ожидание того, как должна выглядеть тень черной дыры (ОТО предсказывает асимметричное кольцо такого-то размера, симуляции GRMHD дают конкретный профиль яркости).

Что произойдет, если крутить ручки регуляризаторов до тех пор, пока картинка не «станет похожа на ожидаемое»? Правильно — она станет похожа на ожидаемое. И в науке такой подход очень рискован.

Команда EHT поступила красиво. Они разбились на четыре независимые группы. Каждая получила одни и те же данные, но не имела права обсуждать промежуточные результаты с другими группами. Группы использовали разные пайплайны и разные настройки гиперпараметров:

• DIFMAP с CLEAN — классика, написанная еще в 90-е;

• eht-imaging — RML на Python, разработанный командой Гарварда;

• SMILI — независимая RML-имплементация, в основном японская;

• Внутри RML-команд — разные комбинации регуляризаторов.

Четыре группы работали параллельно семь недель. Только после этого они встретились и сравнили результаты.

Все четыре получили оранжевый бублик. С одним и тем же диаметром (~40 микросекунд дуги), с одной и той же асимметрией яркости (южная сторона ярче — следствие релятивистского эффекта вращающейся плазмы), с теневой структурой в центре. Это и есть основание утверждать: то, что вы видите на «фото», — не артефакт алгоритма. Четыре независимые имплементации дают один и тот же результат.

Если бы у одной команды получилось одно, у другой — другое, у третьей — третье, то статью бы просто не выпустили.

3,5 петабайт самолетом

Теперь о приземленном. Каждая станция EHT пишет четыре канала (две поляризации × две боковые полосы) с двухбитной оцифровкой. Грубая прикидка потока с одной станции:

4 ГГц × 2 (Nyquist) × 2 бита × 2 поляризации = 32 Гбит/с

С полным стеком из четырех Mark 6 в тандеме станция выходит на 64 Гбит/с агрегатной записи. За семь ночей наблюдений в апреле 2017 года с восьми телескопов так и набрались те самые ~3,5 петабайта.

Передать это по интернету тогда было физически невозможно. Гигабитный канал тащил бы эти данные больше года. Даже 10-гигабитный — больше месяца, а у Южного полюса никаких 10 Гбит и в помине не было (там и сейчас спутниковый узкий канал).

Поэтому данные писали на массивы Mark 6 — стойки с десятками жестких дисков. Стойки физически грузили в самолеты и везли в два корреляционных центра: Хейстекскую обсерваторию (MIT) и Институт радиоастрономии Макса Планка в Бонне.

С Южного полюса диски лететь не могли восемь месяцев — там антарктической зимой ничего не садится и не взлетает. Ученые сидели в Бостоне и Бонне с 7/8 данных и ждали лета в Антарктиде.

Когда диски наконец прилетели, их прогнали через корреляторы — специализированные компьютеры, которые сводят сигналы со всех восьми станций и для каждой пары вычисляют ту самую видимость V(u, v) в куче временных окон. После корреляции 3,5 петабайта сжимаются до условных «гигабайт точек на (u, v)-плоскости с погрешностями». Уже эти гигабайты стали входом для CLEAN и RML.

Кстати, это классический пример теоремы Таненбаума: «никогда не недооценивайте пропускную способность универсала, набитого магнитными лентами и едущего по шоссе». 3,5 петабайта за неделю на самолете — это в среднем около 46 Гбит/с непрерывного потока. Никакой кабель из Антарктиды этого не даст.

Это не фотография

Теперь можно вернуться к началу.

Знаменитое изображение M87 — это численное решение обратной задачи I_dirty = I_true * B_dirty. С регуляризацией, кодирующей физически разумные предположения о небе. На входе у которого — 28 наборов точек на (u, v)-плоскости, полученных корреляцией сигналов с восьми разнесенных антенн. И с предобработкой длиной в два года и пятью петабайтами исходных данных.

Если бы у вас был глаз диаметром с Землю, чувствительный на 1,3 мм, и вы посмотрели бы в сторону Девы — вы увидели бы примерно то же самое. Такого глаза не существует, но его сумели построить, точнее его математический аналог из восьми кусков, разделенных тысячами километров, и решением обратной задачи.

Это не делает картинку менее реальной, просто это другой способ видеть. Половина того, что мы знаем о Вселенной — от первых снимков пульсаров в 60-х до карт реликтового излучения с Planck — получена решением обратных задач разной степени жесткости.

На мой взгляд, EHT — это самый красивый и самый недооцененный пример ,когда четыре независимые команды с четырьмя разными алгоритмами сходятся к одному и тому же результату — это значит, что бублик правда там, в 55 млн световых лет, и весит шесть с половиной миллиардов Солнц.

Вот это и есть фото того года.